ABSTRACT 



VLADIMIROV, ANDREY. Modeling Magnetic Field Amplification in Nonlinear Diffusive 
Shock Acceleration. (Under the direction of Dr. Donald C. Ellison.) 

This research was motivated by the recent observations indicating very strong 
magnetic fields at some supernova remnant shocks, which suggests in-situ generation of 
magnetic turbulence. The dissertation presents a numerical model of collisionless shocks 
with strong amplification of stochastic magnetic fields, self-consistently coupled to efficient 
shock acceleration of charged particles. Based on a Monte Carlo simulation of particle 
transport and acceleration in nonlinear shocks, the model describes magnetic field ampli- 
fication using the state-of-the-art analytic models of instabilities in magnetized plasmas 
in the presence of non-thermal particle streaming. The results help one understand the 
complex nonlinear connections between the thermal plasma, the accelerated particles and 
the stochastic magnetic fields in strong collisionless shocks. Also, predictions regarding the 
efficiency of particle acceleration and magnetic field amplification, the impact of magnetic 
field amplification on the maximum energy of accelerated particles, and the compression 
and heating of the thermal plasma by the shocks are presented. Particle distribution func- 
tions and turbulence spectra derived with this model can be used to calculate the emission 
of observable nonthermal radiation. 
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Chapter 1 

Introduction. Interstellar Shocks, 
Cosmic Rays and Magnetic Fields 

What happens after a massive star explodes at the end of its hfe cycle as a su- 
pernova (SN)? Why are the rims of supernova remnants (SNRs) so thin and luminous in 
the radio, X-ray and gamma ray spectral ranges? Where and how are cosmic rays (CRs) 
produced? What does it take to explain the dynamics of matter in the most energetic sys- 
tems in space, including the cosmological large scale structure of the Universe? The current 
state of affairs in astrophysics makes it clear that, in order to answer these questions, the 
phenomenon of shocks must be studied in detail. The low gas densities in many cosmic 
environments make the shocks collisionless (see Section [l.3p . which gives them properties 
different from those of the collisional terrestrial shocks. 

Understanding shocks is as important for astrophysicists as describing electromag- 
netic waves is for radio engineers. Shocks are born whenever gases or fluids are forced to 
move at a supersonic speed. They compress and heat the interstellar matter (ISM), trans- 
fer energy and momentum, produce cosmic rays that fill and affect the Universe, and, as 
recent observations show, shocks may produce and strongly amplify turbulent magnetic 
fields. Electromagnetic radiation from processes in shocks is a powerful diagnostic of the 
conditions in the shock-generating systems. 
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1.1 Shocks in hydrodynamics 

I would like to illustrate shocks with a phenomenon that we encounter on a daily 
basis - a standing shell shock in a kitchen sink formed by the quickly running water from 
the tap. 




Figure 1.1: Tap water in a sink forms a shock. 

As seen in Figure II. H the falling water hits the bottom of the sink and moves 
outward at a speed that exceeds the speed of the surface waves in water of the local depth. 
This makes a shock form, a relatively stationary enclosed boundary, at which the speed 
of the water flow abruptly drops, and the depth increases. The direction in of the shock's 
apparent motion depends on the choice of the observer's reference frame, but let us adopt a 
convention that unambiguously determines the direction of shock propagation. I will define 
the latter as the direction in which the boundary between the unshocked and shocked 
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media moves with respect to the unshocked medium. In the case of the shock in a sink, the 
unperturbed medium is inside of the circular shell, and it moves outward. Therefore, the 
shock is directed inward (i.e., any small arc of the shock boundary is moving towards the 
center with respect to the water inside the boundary). The arrows show the velocity of the 
water with respect to the shock. 

A similar inward-directed shock exists in the Solar System: the Solar wind, com- 
posed of fast charged particles emitted by the Sun, moves radially outward and collides with 
the cold interstellar material approximately 80-100 AU from the Sun (an astronomical unit, 
1 AU~ 1.5 • 10^^ cm, is close to the distance between the Sun and the Earth). The so-called 
termination shock forms there. At this thin boundary, the Solar wind becomes compressed 
and heated, and its speed drops by a factor of 2-5. Both Voyager spacecrafts recently passed 
through the termination shock on their way out of the Solar System |114^ [24} 1113} 123] . 

1.2 Forward shock of SNRs 

After a star with an initial mass greater than approximately 8 Mq (Mq ~ 2 • 10^^ g 
is the mass of the Sun) runs out of its fusion fuel, or a white dwarf accreting mass from 
another star in a binary system reaches the critical mass and ignites, an explosion will 
occur. This explosion, powered either by gravity, or by thermonuclear fusion, is known as 
a supernova, and ejecting up to 10^^ ergs in kinetic energy, it can be bright enough to see 
with the naked eye thousands of light years away. A remnant of a supernova in our Galaxy 
may remain visible to radio, optical and X-ray telescopes for hundreds or thousands of years 
after the explosion, as it expands into the interstellar medium, cools and gradually fades. 

Hydrodynamic simulations and observations show a common structure of flow that 
forms in SNRs, as shown in Figure [L2l The metal-rich material (ejecta) is thrown out from 
the star at speeds of several thousand kilometers per second. It ploughs through the low- 
density ISM, and eventually forms a strong forward shock in front of it, directed outward. 
A contact discontinuity separates the metal-rich ejecta material from the low-metallicity 
shocked ISM. Simulations show that a reverse shock, may form in the ejecta. While the 
inverse shock is directed inward (i.e., it shocks the material coming from the interior of 
the reverse shock boundary), it may be physically moving outward or inward at different 
stages of the SNR evolution (e.g., [46]). Note that in Figure [L2t the solid arrows show the 
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Figure 1.2: Schematic structure of an SNR. 

velocities of the unshocked medium with respect to the forward and the reverse shocks. The 
dotted hnes indicate the expansion of the forward shock in time. 

Of particular importance to us is the forward shock, because it can be very strong, 
sonic Mach number reaching the values of several hundred. There are two important differ- 
ences between the shell shocks in Figure fTTT] and Figure [TT^ First, the shock in the sink is 
directed inward (it sweeps up water coming from the interior of the circle), while the SNR 
forward shock is directed outward (sweeping up the interstellar matter outside of it). The 
second difference is that the sink shock is stationary, i. e., its radius remains constant in 
time, but the SNR expands into a stationary unshocked ISM, increasing the radius of the 
forward shock. 

One of the hundreds of known and carefully observed SNRs, G 1.9+0.3, stands out 



Figure 1.3: The youngest known galactic SNR, G 1.9+0.3. 



as the youngest known supernova remnant in the Galaxy. It was very recently identified 

as one by an international team led by NCSU astronomers |lU2j and [H7]- It provides an 

illustration of a typical spatially resolved SNR imaged in X-rays. In Figure 11.31 (image 

credit: Prof. S. Reynolds, NCSU, [102J), the dotted line maps the approximate location of 

the forward shoclo; the dotted arrows indicate the direction of the shock movement with 

respect to the interstellar medium, and the solid arrows indicate the directions and the 

relative magnitudes of the velocities of the unshocked (the long arrow) and the shocked 

(the short arrow) plasma, with respect to the shock. Note that in the following text we 

usually adopt the reference frame in which the shock is at rest, and the plasma is flowing 

^The determination of the location of the shock is a comphcated problem, and the contour in Figure [L3l 
should be perceived as an artist's impression. 
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into the shock at a supersonic speed. This approach corresponds to the flow directions 
shown in Figures 11.11 and II. 3[ 

1.3 The concept of a collisionless shock 

Generahy, gas flowing into a shock gets compressed and heated in a narrow region. 
But how narrow can this region be for a shock in an astrophysical plasma? In order to change 
the density, bulk speed and temperature, the gas particles must experience a few strong 
collisions, and the thickness of the shock can therefore be estimated as the mean free path 
of particles between collisions. Indeed, for shocks in dense gases (for example, air) a particle 
mean free path is comparable to the shock thickness. However, in an attempt to apply the 
same reasoning to interstellar or interplanetary shocks, one runs into a complication. 

For a plasma consisting of fully ionized hydrogen, the cross section of Coulomb 
collisions between protons is formally inflnite [110], but we can roughly estimate the cross 
section of collisions that are strong enough to change the energy of the particles significantly. 
If by 'significantly' one means that the change energy due to collision must be comparable 
to the thermal energy, then the protons must approach each other within a distance Tc such 
that 

- = ksT. (1.1) 

Here and in the rest of the equations in this dissertation, the CGS system of units is 
adopted. The quantity e is the elementary charge, and the left-hand side of Equation (jl.ip 
is the electrostatic potential energy of two protons separated by the distance rc- The right- 
hand side is the characteristic thermal energy of protons in a gas of temperature T (ks is 
the Boltzmann constant). This gives a rough estimate of the collision cross section 

4 

^ = = (1-2) 



and of the mean free path 

A=- = ^, (1.3) 



an vre^n 



where T is the temperature and n is the number density of the gas. For conditions typical 
for the Solar Wind in the near Earth space, n ~ 4 cm~^, T ~ 10^ K, which gives A ~ 
3 • 10^^ cm= 2 • 10^ AU. This distance is much greater than the size of the Solar System, 
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which means that shocks just do not have room to form in the Solar wind near the Earth. 
However, spacecraft observations clearly indicate numerous interplanetary shocks of various 
strengths traversing the Solar System. Measurements reveal that the interplanetary shocks 
are much thinner then the number above: observed thicknesses are around A ~ 10^ — 10^'^ cm 

m- 

These observational data are successfully explained by the theory of collisionless 
shocks, which assumes that in the transition region of the shock, particles collide not with 
each other, but with inhomogeneities of magnetic fields. This shrinks the thickness of the 
transition region down to the scales of (multiple) proton gyroradii (see Section 6.4 of [81j). A 
shock in which collisions between particles play a negligible role compared to the dynamics 
of the particles in stochastic magnetic fields is called a collisionless shock, and the term 
collisionless plasma is widely used to define the systems in which similar conditions exist. 

An interesting property of collisionless plasmas is that, due to the absence of 
particle-particle collisions, the time scales of thermalization of non-equilibrium energy dis- 
tributions of particles are extremely large. This allows for the existence and sustainability 
of a superthermal component in the particle distribution (i.e., energetic particles). Present 
research, along with other models, shows that the superthermal particles may be not just a 
minor admixture to the thermal particle pool, but, on the contrary, they may dominate the 
dynamics of a collisionless shock. This assertion is explained in the following two sections. 

1.4 Cosmic rays 

Cosmic rays (CRs) are charged particles, first seen as radiation coming from space 
in a balloon experiment performed by Victor Hess in 1912, and identified as charged nuclei by 
Phyllis Freier and others in 1948 |57]- The spectrum of these particles spans many decades 
in particle energy (lO"^ to lO^'' eV (!) per nucleus) as well as in flux (from 1 cm^^s^^sr^^ 
for energies of 1 GeV and above, down to 1 particle per square kilometer per century for 
energies over 10^^ eV j34|). 

From the multitude of observational data on CRs, it is known that the lower energy 
CRs come from the Sun, and the higher energy CRs (over 1-10 GeV) are of Galactic origin. 
CRs are therefore the second most important source of information about deep space after 
electromagnetic radiation. The problem of measuring and explaining the spectrum, compo- 
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sition, temporal variation and directional distribution of CRs is extensive and longstanding. 
It requires answering two major questions: how CRs are produced, and what happens to 
them en route from the source to the detector on Earth. 

The spectrum shown in Figure [ITS] (image credit: S. Swordy, University of Chicago, 
[116| ) is a compilation of the various measurements. This spectrum cannot be identified 
with a single CR source or even multiple CR sources; in fact, it represents a superposition 
of the multitude of Galactic CR sources, integrated over a time scale of millions of years, 
and convolved with the history of their propagation in the Galactic magnetic fields from all 
parts of the Galaxy. 

Most researchers these days are convinced that the bulk of the Galactic CRs at 
least up to the 'knee' of the CR spectrum (i.e., up to the energies of 3-10^^ eV), are produced 
in astrophysical collisionless shocks p3]. While the details of this process may be uncertain 
(how much shocks of individual SNRs contribute to the CR production, in comparison with 
SNR shock ensembles in the so-called superbubbles, how the interstellar dust is involved, 
etc.), the general idea is commonly accepted now. 

The process that accelerates the particles to ultra-relativistic energies in shocks is 
known as the first order Fermi procesj^, or diffusive shock acceleration (DSA). 



1.5 DSA — test-particle approximation 

Diffusive shock acceleration, DSA, also known as the first order Fermi process 
(often abbreviated as Fermi-I) was first applied to the problem of cosmic ray production in 
shocks by several independent groups and researchers at the end of the 1970s [9} [20l [79} [30] . 
The best simple mechanical analogy to this process is the acceleration of a rubber ball 
elastically bouncing back and forth between two massive walls, as the walls are slowly 
moved towards each other. In a shock, the role of the moving walls is played by the bulk 
gas flow: the faster-moving unshocked gas and the slower moving shocked gas form an 
effectively converging system. 

The spectrum of particles accelerated in such a manner may be calculated in 
various ways, including a kinetic approach (see, e.g., [TOj). Consider a one-dimensional 

^There also exists a model of the second order Fermi process, in which particles are accelerated by 
magnetohydrodynamic turbulence in the absence of shocks. 
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Figure 1.4: The all particle spectrum of cosmic rays 
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shocked flow, with the shock located at x = 0, and the flow speed 

u{x) = { ^ (1.4) 




where uq is the upstream and U2 < uq - the downstream speed, and let there be a minor 
admixture of energetic particles that move diffusively in the bulk plasma, the diffusion 
being isotropic in the plasma frame. Assume that the diffusion coefficient is independent of 
momentum and of coordinate (except it may have different constant values upstream and 
downstream of the shock): 

-Do, X <0, 



D{x) 



D2, X > 0, 



In a steady state, diffusive propagation of the energetic particles, as they are being advected 
downstream by the flow, can be described by the equation 

„w^a|^=p(x)%i2, (1.6, 

where f{x,p) is the particle distribution function, such that f (x , p)dxdydzdpxdpydpz is the 
number of particles in the phase space volume dxdydzdpxdpydpz, and that f{x,p) docs not 
depend on the direction of p. Suppose the incoming energetic particles have a distribution 
function: 

lim f{x,p) = fo{p) = fo^Soip - Po), (1.7) 

X-+-00 p^ 

where po is a momentum such that the corresponding particle speed is much greater than 
no, p is the current particle momentum, and Sd is the Dirac delta- function. Assume the 
trivial downstream boundary condition 

lim f{x,p) < 00, (1.8) 

and define the conditions at the discontinuity point x = 0: 

lim f{x,p) = lim f{x,p), (1.9) 

(_D^df^_p_df{^\ ^ lin. (_D^df{x^P)_Pjny))_ (i.ip) 
x-*o-\ " 5x 3 dp J X-.0+V ^ dx ?, dp J ^ ' 
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The first equation expresses the requirement of continuity of the particle density, and the 
second ~ of particle flux. The general solution of equation (jl.6p may be written as 



fix,p) = < 



C{p)exp(^j+E{p), x>0. 



(1.11) 



Substitution of this form into the boundary condition (|1.7p results in 

Bip) = hip), 



[1.12) 



and using the boundary condition (jl.Sp gives 



(1.13) 



C{p) = 0. 

Now we can use the conditions at x = 0, where the density continuity equation (jl.Op can 
help constrain A[p) and E{p) in (|l.lip : 



A{p) + Up)=E{p), 
and flux continuity condition (jl.lOj) . rewritten as 

0„ l»J^\-D, lim (m^\^_v(m>:A 



dx 



' x^o+ V dx 



3 V dp 



gives 



DoA{p)^ exp (0) - = -PJEM. _ ^2) , 
Uq 3 dp 



Combining (jl.l4p and (jl.l6p . we get 

p dE{p) 



fo 



{uq-U2) + E{p)ui=ui^5d{p-Pq), 
6 up Pq 

which can easily be integrated, assuming £'(0) = 0, and the solution is 

Eip)=Eo(^^y H{p-po), 



where 



En 



Pq [uo - U2) 
3tto 

Mo - ""2 ' 



(1.14) 

(1.15) 
(1.16) 

(1.17) 

(1.18) 

(1.19) 
(1.20) 
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and H{z) is the Heaviside step function: 

H{z) = { " - ^ (1.21) 




Finally, the solution of equation (II. 6p with boundary conditions (II. 7p . (jl.Sp and the conti- 
nuity conditions at the shock (jl.Op and (jl.Op is: 



VP/ Po ^ ' (1.22) 

i^ofyj Hip-po), x>0. 

This is the so-called test-particle solution of the problem of diffusive shock acceleration, 
meaning that the accelerated energetic particles are implicitly assumed to be a small ad- 
mixture in the vast thermal pool. This assumption is likely to fail for strong collisionless 
shocks, leading to serious modifications of the solution, on which the present work concen- 
trates. 

Let us analyze the basic properties of the test-particle solution (jl.22p . 

• It requires that some seed particles be introduced, represented by fo{p), but in real 
shocks these seed particles must be produced from the thermal pool {injected, as the 
theorists of the particle acceleration field prefer to put it). This model is unable to 
predict anything about the injection of particles, and their number /o and momentum 
Po are free parameters of the test-particle model. 

• Once the seed particles are introduced, they form a power-law superthermal tail up- 
ward of the injection momentum po, with the index s that depends only on the pre- 
shock and the post-shock speed, as given by equation 11.201 That equation can be 
re- written in terms of the shock compression ratio r = uq/u2 as 

s = = -. (1.23) 

UQ-U2 r - 1 

For the strongest hydrodynamic shocks in a non-relativistic monatomic gas, the com- 
pression ratio r = uq/u2 approaches the value of r = 4 (this well known result can 
easily be derived from the Hugoniot adiabat presented in Section 13.1.71 in the limit 
Ms — > oo with 7 = 5/3) . Notably, such compression ratio corresponds to the power law 
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index of the accelerated particle distribution s = 4. A particle distribution f{p) oc 
extending to p ^ cxd in unphysical, because the internal energy of such distribution 
diverges logarithmically at p — > oo. This means that, if compression ratios of r = 4 
or greatecl are achieved in space, there must be some process responsible for limiting 
the maximum achievable energy. The escape of the highest energy particles from the 
system, or a finite time of particle acceleration in a time-dependent calculation may 
determine the high-energy cutoff of the particle spectrum. Such processes are not 
included in this simplistic model. 

• The basic physical assumption that leads to the emergence of the power-law superther- 
mal tail of f{p) is that the particles are subject to diffusion isotropic in the plasma 
frame (this is expressed by the equation (II. 6p ). This implies that we are dealing with 
a collisionless shock (otherwise the superthermal particles would have to thermal- 
ize through collisions with their thermal counterparts) that has a certain stochastic 
magnetic field structure (i.e., turbulence), responsible for particle scattering. The 
properties of these stochastic fields are, obviously, beyond the scope of this model, 
but they must influence the solution by at least determining the diffusion coefficient 
D{x,p). In fact, as will be shown later, the magnetic turbulence that confines the 
particles to the acceleration site, and allows for the Fermi-I acceleration, is probably 
produced by the accelerated particles themselves, which raises the question of solving 
the particle acceleration problem consistently with the turbulence production process. 

It turns out that all of these properties of the test-particle make it unable to explain some 
observations of interstellar collisionless shocks (see the next section) , which calls for a better 
model. 



^For example, relativistic gases with a polytropic index 7 = 4/3 allow the strongest shocks to have a 
compression ratio up to r = 7, which results in a power law index s = 3.5 - an even more strongly diverging 
distribution. 
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1.6 DSA — nonlinear regime 

The example from classical mechanics that illustrates the first order Fermi process 
- the rubber ball bouncing between two converging walls - may also be used to understand 
the nonlinear aspects of diffusive shock acceleration. The ball in classical mechanics gains 
energy at every collision, but only as long as the walls are much heavier than the ball and 
continue to move inward despite the ball's kicks. But what if the ball gains enough energy, 
so the recoil of the walls makes them slow down their convergence? In this case we would 
have to account for the feedback of the ball on the walls. This makes the problem nonlinear. 
Suppose, we put one ball between the walls, and its kinetic energy after cycles becomes 
K. If we were to put not one, but two balls, the their total energy after N cycles would be 
less than 2K, because the recoil of the two balls would have slowed the walls down more 
efficiently than the recoil of one ball. Similarly, in shock acceleration, the energy of the bulk 
plasma flow powers the energetic particle acceleration, but once the accelerated particles 
gain enough energy to push back on the flow, the situation changes dramatically. Such a 
system is called a nonlinear, or a multicomponent shock wave, and is the subject of study 
of nonlinear diffusive shock acceleration theory. 

There is a number of reasons to believe that strong shocks in space accelerate 
particles very efficiently, thus operating in the nonlinear regime. I outline these reasons 
below, and some of them will be elaborated on further in the dissertation. 

1. Energy considerations. The energy density of Galactic cosmic rays at the location of 
the Solar System is £ 0.6 eV-cm~^, and their characteristic age inferred from the 
radioactive nuclei in CRs is of order Tcr ~ 10^ yr. Assuming that the escape of CRs 
from the Galaxy has the time scale Tcr and that the escape is balanced by the CR 
production one can estimate the required power of CR production in the Galaxy as 

eV 

Per = — = 3-10^^ erg s-\ 

where V is the volume of the Galaxy, V = 7r(10^ pc)^ • 100 pc = 3 • 10^^ pc^ ^ 10^^ cm^. 
Assuming that all of these cosmic rays are produced by shocks of SNRs, which occur 
once in Tgn = 100 yrs and release E = 10^^ erg as the kinetic energy of the shock 
wave, the Galactic energy production in the form of shocks is 

Psk = — = 3-10^^ erg s-\ 
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Our estimates of the quantities Per and -Psk are comparable, which means that SN 
shocks may easily be required to have an efficiency on the order of tens of percent of 
converting the bulk motion energy into the energy of accelerated particles. See also 
|19j for a detailed discussion. 

2. Numerical simulations (e.g., [l3l[89l[711[22]) predict efficient particle acceleration given 
the simplest physically realistic model of particle injection, thermal leakage. The above 
mentioned models use different techniques, but all of them predict that in a strong 
collisionless shock, the energy density of energetic particles becomes comparable to 
the kinetic energy density of the flow, thus making the problem nonlinear. 

3. Analysis of the morphology of resolved SNRs indicates high compression ratios at the 
forward shock (e.g., [125[ I33j ) which is consistent with the predictions of nonlinear 
particle acceleration theories. There is also observational evidence of the nonlinearity 
of DSA that comes from the analysis of nonthermal emission spectra (e.g., [2]) 
and from spacecraft observations of the Earth's bow shock (e.g., [51j). 

4. Recent observations indicate that magnetic fields in some SNR shocks are much 
stronger than the ambient magnetic fields, which makes many researchers believe 
that magnetic fields are amplified in situ, i.e. in the shocks, by the shock-accelerated 
particles. If that is the case, the energy density of the accelerated particles must be 
no less than the energy density of the amplified magnetic fields, and, according to the 
observational estimates, the latter is a significant fraction of the dynamical pressure 
of the shock flow. This necessitates the nonlinear DSA (see, e.g., [53])- 

The nonlinear DSA theory was developed by various researchers in the 1980s. 
Although the details of the models may differ, all of them agree on the following: when 
the energetic particles gain enough energy to feed back on the flow, the unshocked plasma 
slows down and becomes compressed even before it reaches the viscous shock (the latter 
is renamed a subshock in context of nonlinear DSA), which means that a shock precursor 
forms in the upstream region (x < 0); the maximum particle energy must be limited either 
by the age, or by the size of the shock, and if particles of the highest energies are allowed 
to leave the shock far upstream, it leads to an increase in the compression ratio. 
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1.7 Magnetic field amplification in shocks 

Recent observations and modeling of several young supernova remnants (SNRs) 
suggest the presence of magnetic fields at the forward shock (i.e., the outer blast wave) well 
in excess of what is expected from simple compression of the ambient circumstellar field, 
i?ism- These large fields are inferred from: 

• spectral curvature in radio emission (e.g., [IB]) , 

• broad-band fits of synchrotron emission between radio and non-thermal X-rays (e.g., 
[IllISl], see also l35j), 

• sharp X-ray edges (e.g., [Ml H [IMl El [32] ) , and 

• rapid variability of nonthermal X-ray emission from bright filaments in SNRs (first 
reported by [119] ). 

While these methods are all indirect, fields greater than 500 /iG are inferred in the supernova 
remnant Cassiopeia A and values of at least several 100 fiG are estimated in Tycho, Kepler, 
SN1006, and G347.3-0.5. If -Bism ~ 3 — 10 /uG, amplification factors of 100 or more may 
be required to explain the fields immediately behind the forward shocks and this is likely 
the result of a nonlinear amplification process associated with the efficient acceleration 
of cosmic-ray ions via diffusive shock acceleration (DSA). The magnetic field strength is 
a critical parameter in DSA and also strongly infiuences the synchrotron emission from 
shock accelerated electrons. Since shocks are expected to accelerate particles in diverse 
astrophysical environments and synchrotron emission is often an important emission process 
(e.g., radio jets), quantifying the magnetic field amplification has become an important 
problem in particle astrophysics and has relevance beyond cosmic-ray production in SNRs. 

These highly amplified magnetic fields are most likely an intrinsic part of efficient 
particle acceleration by shocks. This strong turbulence, which may result from cosmic ray 
driven instabilities, both resonant and non-resonant, in the shock precursor, is certain to 
play a critical role in self-consistent, nonlinear models of strong, cosmic ray modified shocks. 
Although plasma wave instabilities in presence of accelerated particles have been studied in 
the context of shock acceleration before (e.g., [HI [82]), it was only recently suggested that 
these instabilities may lead to very efficient amplification of magnetic field fiuctuations. 
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AB S> Bq \TT\. Since then, new models of plasma instabilities possibly responsible for 
efficient magnetic field amplification were proposed [861 lOl EO] and studied in context of 
shock acceleration [991131 [T22l [T29] . 

All these plasma instabilities are assumed to amplify pre-existing waves in a plasma 
in the presence of an underlying uniform magnetic field Bq parallel to the flo'VMj. The two 
models of interest that will be applied to the present work are: 

1. Resonant CR streaming instability (see [106^ [82] and [TTl [3l [T22] ). in which particles 
of a certain momentum amplify Alfven waves with a wavenumber equal to the inverse 
gyroradius of the particle, and 

2. Nonresonant CR streaming instability of short- wavelength modes (suggested by [lOj). 
which I will sometimes refer to as Bell's instability, in which the diffusive electric 
current of CRs amplifies almost purely growing waves with wavenumbers much greater 
than the inverse particle gyroradius. 

I should also mention a nonresonant instability that produces long-wavelengths modes and 
may also be important in shocks (see [30j and Section [3. 2. 3p . which I am planning to apply 
to the modeling of shocks in the future, as well as other possible mechanisms (e.g., [9l]). 

Amplification of magnetic turbulence has great importance in the process of DSA. 
The amplified turbulence provides the stochastic magnetic fields that scatter the acceler- 
ated particles, allowing them to participate in the Fermi-I process. The properties of the 
particle scattering are therefore dependent on the spectrum of stochastic magnetic fields, 
yet the latter are produced by the accelerated particles. This complex connection between 
particles and waves in shocks adds to the nonlinear nature of shock acceleration, discussed 
in Chapter [2l Therefore, magnetic field amplification affects the observable nonthermal 
synchrotron emission from shocks in two ways: it determines the structure and strength of 
the magnetic fields in which the emission occurs, and shapes the spectrum of the radiating 
energetic particles. 



*For relativistic shocks, which are beyond the scope of the present research, the Weibel instabihty may be 
important for magnetic field generation, as long as the upstream magnetic field is low (see, e.g., [95II87I[109] 1. 
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1.8 Turbulence 

Studying strong magnetic field amplification in interstellar shocks inevitably makes 
us face the subject of magnetohydrodynamic (MHD) turbulence. Usually turbulence is 
defined as chaotic fluid motion, that is, a motion with a very sensitive dependence on initial 
conditions. Chaotic behavior makes turbulent motions effectively non-deterministic, but 
they can be studied using statistical methods. 

Motions of gases and fluids of high Reynolds number tend to transit to the tur- 
bulent regime (see, e.g., \83\ 196]). which is encountered on a regular basis in areas ranging 
from plasma fusion engineering and race car design to air transport, plumbing, golf and 
food processing (e.g., [581 IllSj [TBI I112j ). Driven by the need of applications like meteo- 
rology, climate modeling, aerospace engineering, and others, turbulence research has been 
conducted for many decades, and is a challenging field of mathematics and physics (e.g., 
[56)). Conducting fluids (plasmas) easily develop and sustain magnetic fields, and the MHD 
turbulence regime, occurring in plasmas, is even more complicated by the magnetic field 
interactions than its hydrodynamic counterpart |17j . 

Considering that plasmas constitute a large fraction of all baryonic matter in 
space, their properties have pervading importance for astrophysics. Namely, turbulence 
in plasmas determines cosmic ray acceleration and propagation, plays a crucial role for 
angular momentum transfer in accreting systems and impacts the properties of gravitational 
collapse. The list of astrophysical objects affected by MHD turbulence is therefore extensive: 
large scale structure of the Universe, quasars, accreting binary systems, forming stars, 
supernova remnants, etc. 

The primary sources of information about MHD turbulence are spacecraft obser- 
vations of interplanetary space and numerical simulations. The former provide real, but 
often hard to interpret data, the bottom line of which is that turbulence often consists of 
stochastic perturbations of plasma velocities and magnetic fields spanning many decades 
of the spatial scales. Oftentimes, the Fourier spectrum of spatial structure of turbulent 
fluctuations reveals a power-law distribution of energy in wavenumber space. The numeri- 
cal simulations have the advantage of providing data that is easy to analyze and scale for 
practically applicable theories. 

A simplified picture of turbulence evolution, based on extensive research, involves 



19 



three dominant processes: energy supply, spectral transfer of energy and dissipation. Con- 
sider a fluid flow in a pipe, where a large flux of the fluid leads to the development of a 
hydrodynamic instability that creates vortices (eddies) breaking the laminar flow. In this 
way energy is supplied to the turbulence in the form of large-scale vortices. These eddies 
then break down into smaller eddies ~ this way, spectral energy transfer (cascade) from 
large to small scales is realized. As the scale of the turbulent structures due to cascading 
becomes smaller, fluid viscosity plays an increasingly greater role, eventually leading to the 
dissipation of the smallest eddies into heat. 

The MHD turbulence, as mentioned above, is difficult to describe. It was originally 
treated and analyzed as a set of small perturbations (i.e., plasma waves) moving in the 
large-scale uniform magnetic field and weakly interacting with each other (the so-called 
Iroshnikov-Kraichnan approach |70|, I78j). However, Goldreich and Sridhar [64 



^' point out 



that this approach may be inappropriate for MHD turbulence due to its inherent anisotropy 
introduced by the magnetic field [55]. The bottom line of their theory and of the subsequent 
simulations of MHD is that the magnetic field plays a stabilizing role. The cascading takes 
place mostly for wave vectors perpendicular to the uniform magnetic field, while the parallel 
cascade is suppressed. 

A comprehensive source on classical theory of hydrodynamical turbulence is [96] . 
Modern advances in the study of MHD turbulence is presented in [T7] . 



^Note that this publication is relatively recent, but it has laid the groundwork for MHD turbulence 
research from a new vantage point, possible only with modern computational resources. 
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Chapter 2 

The Problem of Nonlinear DSA 

The general problem of nonlinear diffusive shock acceleration of charged particles 
(DSA) can be formulated as follows: given a supersonic flow with a speed uq of a plasma 
with a number density no, temperature Tq and a pre-existing magnetic field Bq, and given 
the location x = where this flow develops a subshock, find the distribution of particles 
/(x, p, t) and electromagnetic fields in the shock vicinity. This problem is complicated 
by two facts: a) particle acceleration occurs due to complex motions of particles in the 
turbulent magnetic field, but the magnetic turbulence itself is dependent upon the motion 
of the accelerated particles, and, b) if particle acceleration is efficient, different parts of the 
particle spectrum interact with each other (i.e., the accelerated particles push back on and 
slow down the flow of the thermal particles). 

This problem cannot be practically tackled by particle simulations from first prin- 
ciples like Maxwell's equations and Lorentz force (see Section [2. 2p . and the most computa- 
tionally expensive operations must be performed analytically. Namely, all currently existing 
models of nonlinear DSA, including the one discussed in this dissertation, assume that the 
accelerated particles propagate diffusively with some diffusion coefficient, or mean free path, 
prescription. This allows the models to eliminate the need to describe the complex inter- 
actions between particles and waves, and to concentrate on the physical aspects of particle 
acceleration. 
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2.1 Analytic models 

In successful analytic models, a one-dimensional steady state shock with a non- 
linear precursor is described by the flow speed u{x,t), mass density of the plasma p{x,t) 
and an isotropic distribution function of energetic particles, fcT{x,p,t). The above men- 
tioned macroscopic quantities must be consistent with the fundamental conservation laws: 
mass, momentum and energy must be conserved. These conditions are expressed with the 
following system of equations: 

pu = const, (2.1) 

/OU^ + Pth + -Per + -Pmag = COnst, (2.2) 
\pu^ ^ ^^^^-Pth^tH -^^^-Per^i + ^-Pmaglt + Qesc = COnst. (2.3) 

2 7th - 1 7cr - 1 2 



and the evolution of the particle distribution is governed by the kinetic equation of CR 
transport 



d_ 

dx 



d 

D{x,p)—f{x,p) 



, d f{x,p) _^ 1 / ^\ df{x,p ) 
dx (i \dx J dp 



Equations (j2.ip . (j2.2p and (j2.3p represent conservation of mass, momentum and energy 
fluxes, respectively. Equation (j2.4p is the kinetic equation describing propagation of cos- 
mic rays in the diffusion approximation. The expressions above are, essentially, a direct 
generalization of the test particle model of shock acceleration demonstrated in Section 11.51 
complemented by the treatment of the flow speed u{x) variability upstream. Let us use 
the following notation for the flow speed and other quantities at points of interest: uq is 
the far upstream flow speed, U2 is the downstream flow speed, and ui is the flow speed 
just before the subshock. Thus, in the upstream region, x < 0, the flow speed varies from 
u{x = — oo) = uq to u{x = —0) = ui < uo, and then jumps in a viscous subshock to 
u{x = +0) = U2 < ui. Let us also define the total compression ratio, rtot = uo/u2 and the 
subshock compression ratio r^nh = Ui/u2- 

To close the model, one must describe the evolution of thermal gas pressure, Pth) 
define the cosmic ray pressure. Per and have a model for determining the magnetic field 
pressure, -Pmag- The term Qesc representing the energy escape from the system requires 
a model of particle escape, the term Qinj representing the injection of thermal particles 
into the acceleration process calls for a model of particle injection. The most important 
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parameter of the model, the diffusion coefficient D{x,p), must be calculated using some 
simple approximation or using the assumed spectrum of turbulence. 

These models were used in [12l |89l [75l [22l Hi]. The major advantage of the 
analytic models is that they provide a fast solution of the problem, which can be used in 
the simulations of objects incorporating shocks, for example, to calculate the spectrum of 
electromagnetic emission from an evolving SNR (see, e.g., [16^ I97j). 

The computation speed comes at the cost of making some important approxima- 
tions, as summarized below. 

1. Analytic models are limited by the assumptions that go into the analytic description 
of the plasma physics. For example, the diffusion coefficient of charged particles in 
stochastic magnetic fields can be reliably estimated either in the limit of weak turbu- 
lence, or in the simplistic Bohm approximation. Similarly, the analytic description of 
the physics of turbulence generation is only valid in the quasi-linear regime, i.e., weak 
turbulence. 

2. These models adopt a diffusion approximation of particle transport. Hand in hand 
with this approximation goes the assumption that the particle distribution function is 
isotropic in the plasma frame. Only this way can one define the pressures of thermal 
particles Pthi^, t) and cosmic rays -Per (2;, t) as moments of particle distribution function 
f{x,p,t). The isotropy assumption breaks down for relativistic shocks. But even 
for the non-relativistic shocks that we are discussing in this work, the anisotropy of 
particle distribution is important, because it determines the thermal particle injection 
process. Therefore, analytic models require additional parameters or assumptions in 
order to estimate particle injection. 

3. If a strong uniform magnetic field is present in the shock, its strength and orientation 
may affect particle injection and transport. Some SNRs may have an asymmetric 
appearance due to the variation of the obliquity of magnetic field around the rim |33j . 
Analytic models are not able to account for this effect due to the isotropy assumption. 

4. Including some important physical processes in the analytic models of shock acceler- 
ation complicates calculations. These important processes include particle escape far 
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upstream, nonlinear processes in turbulence generation (such as cascading), modifi- 
cations of the diffusion regime by a specific shape of turbulence spectrum, etc. 

Despite these limitations, analytic models are very helpful for making qualitative 
and quantitative predictions regarding the nonlinear structure of shocks, and are the current 
method of choice for modeling the electromagnetic emission of SNRs, where fast simulations 
of shocks are required. 

2.2 Particle-in-cell (PIC) codes 

In principle, the problem can be solved completely with few assumptions and 
approximations with plasma simulations. Those fall into two major categories: particle-in- 
cell (PIC) simulations (e.g., [108^ [98]). and hybrid models that assume that electrons are 
not dynamically important (e.g., [1271 161j j^. 

However, modeling the nonlinear generation of relativistic particles and strong 
magnetic turbulence in collisionless shocks is computationally challenging and PIC simula- 
tions will not be able to fully address this problem in nonrelativistic shocks for some years to 
come even though they can provide critical information on the plasma processes that can be 
obtained in no other way. In this section I outline the requirements that a PIC simulation 
must fulfill in order to tackle the problem of efficient DSA with nonlinear magnetic field 
amplification (MFA) in SNR shocks. The reasoning presented in this section was also laid 
down in |123j . 

There are two basic reasons why the problem of MFA in nonlinear diffusive shock 
acceleration (NL-DSA) is particularly difficult for particle-in-cell (PIC) simulations. The 
first is that PIC simulations must be done fully in three dimensions to properly account for 
cross-field diffusion. As Jones [73j proved from first principles, PIC simulations with one 
or more ignorable dimensions unphysically prevent particles from crossing magnetic field 
lines. In all but strictly parallel shock geometry!^ a condition which never occurs in strong 
turbulence, cross-field scattering is expected to contribute importantly to particle injection 

^ We must also mention the MHD models (e.g., [10] . [130j . [129j l that ignore or treat in a simplified 
way the spectral properties of the particle distribution; while they may be important for describing certain 
aspects of plasma physics, their application to nonlinear DSA is limited 

^Parallel geometry is where the upstream magnetic field is parallel to the shock normal. 
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and must be fully accounted for if injection from the thermal background is to be modeled 
accurately. 

The second reason is that, in nonrelativistic shocks, NL-DSA spans large spatial, 
temporal, and momentum scales. The range of scales is more important than might be 
expected because DSA is intrinsically efficient and nonlinear effects tend to place a large 
fraction of the particle pressure in the highest energy particles. The highest energy particles, 
with the largest diffusion lengths and longest acceleration times, feed back on the injection 
of the lowest energy particles with the shortest scales. The accelerated particles exchange 
their momentum and energy with the incoming thermal plasma through the magnetic fluc- 
tuations coupled to the flow. This results in the flow being decelerated and the plasma 
being heated. The structure of the shock, including the subshock where fresh particles are 
injected, depends critically on the highest energy particles in the system. 

A plasma simulation must resolve the electron skin depth, c/upe, i-e., i^ceii < c/upe, 
where Upe = [47rnee^/me]^/^ is the electron plasma frequency and Lceii is the simulation cell 
size. Here, is the electron number density, rrie is the electron mass and c and e have their 
usual meanings (the speed of light and the elementary charge, respectively) . The simulation 
must also have a time step small compared to ojp^, i.e., ttstep < '^pe- wishes to follow 

the acceleration of protons in DSA to the TeV energies present in SNRs, one must have 
a simulation box that is as large as the upstream diffusion length of the highest energy 
protons, i.e., K(Emax)/^io ~ ''9(-E'max)c/(3uo), where k is the diffusion coefficient, rg{Ejaa.x) 
is the gyroradius of a relativistic proton with the energy Emax^ the shock speed, and 

assuming Bohm diffusion. The simulation must also be able to run for as long as the 
acceleration time of the highest energy protons, Tacc(^max) ~ -E'maxc/(ei?ttQ). Here, B is 
some average magnetic field. The spatial condition gives 

{cj^p,) VTeVy VlOOOkms-V Vcm-3y' ^1836^ ' ^ 

for the number of cells in one dimension. The factor / = mp/rrif, is the proton to electron 
mass ratio. From the acceleration time condition, the required number of time steps is, 

^^-6-10 (^Wj VlOOOkms-0 VT^J -^ ^ ^ 

Even with / = 1 these numbers are obviously far beyond any conceivable computing capa- 
bilities and they show that approximate methods are essential for studying NL-DSA. 
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One approximation that is often used is a hybrid PIC simulation where the elec- 
trons are treated as a background fluid. To get the estimate of the requirements in this case, 
we can take the minimum cell size as the thermal proton gyroradius, rgo = c^2myE<^\^j (eB). 
Now, the number of cells, again in one dimension, is: 

^-'^ IwJUoOOkms-iJ Ikivj • ^'-^^ 

The time step size must be Tstep < ^cp^ where Wcp = eB/nipC is the thermal proton 
gyrofrequency. This gives the number of time steps to reach 1 TeV, 

''"acc(-S'max) ^ ^ f -^max \ / ''^O \ ^2 g~\ 

Ucp ^ ' V TeV y VlOOOkms-V ' 
These combined spatial and temporal requirements, even for the most optimistic case of 
a hybrid simulation with an unrealistically large Tstep, are well beyond existing computing 
capabilities unless a maximum energy well below 1 TeV is used. 

Since the three-dimensional requirement is fundamental and relaxing it eliminates 
cross-field diffusion, restricting the energy range is the best way to make the problem ac- 
cessible to hybrid PIC simulations. However, since producing relativistic particles from 
nonrelativistic ones is an essential part of the NL problem, the energy range must com- 
fortably span mp(? to be realistic. If -Bmax = lOGeV is used, with uq = 5000km s^^, 
and ^th = lOMeV, equation (f27fl) gives ~ 1400 and equation (fZ^Sj) gives ~ 4 • 10"^. Now, 
the computation may be possible, even with the 3-D requirement, but the hybrid simula- 
tion can't fully investigate MFA since electron return currents are not modeled. The exact 
microscopic description of the system is not currently feasible. 

It's hard to make a comparison in run-time between PIC simulations and the 
Monte Carlo technique used here because we are not aware of any published results of 
3-D PIC simulations of nonrelativistic shocks that follow particles from fully nonrelativistic 
to fully relativistic energies. A direct comparison of 1-D hybrid and Monte Carlo codes 
was given in [17] for energies consistent with the acceleration of diffuse ions at the quasi- 
parallel Earth bow shock. Three-dimensional hybrid PIC results for nonrelativistic shocks 
were presented in [62] and these were barely able to show injection and acceleration given 
the computational limits at that time. As for the Monte Carlo technique, a simulation 
that calculates the nonlinear structure of a shock with a dynamic range typical for SNRs, 
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typically takes several hours on 4-10 processors. Thus, realistic Monte Carlo SNR models 
are possible with modest computing resources. 

Despite these limitations, PIC simulations are the only way of self-consistently 
modeling the plasma physics of collisionless shocks. In particular, the injection of thermal 
particles in the large amplitude waves and time varying structure of the subshock can only 
only be determined with PIC simulations (e.g., [98l I108j ). Injection is one of the most 
important aspects of DSA and one where analytic and Monte Carlo techniques have large 
uncertainties. 

2.3 Monte Carlo Simulation 

The Monte Carlo method of solving the problem of nonlinear DSA was devel- 
oped by Ellison and co-workers (see |49t 1721 1122j and references therein for more complete 
details). This method provides an excellent compromise between the physically realistic, 
but computationally limited PIC simulations and fast, but simplified analytic models. The 
compromise is achieved by replacing the solution of coupled equations of particle propa- 
gation, of conservation laws, and of turbulence generation, with a Monte Carlo simulation 
of particle transport that incorporates an iterative procedure that ensures the simultane- 
ous consistency of all assumed laws. The Monte Carlo method goes beyond the diffusion 
approximation of particle transport and allows the calculation of rates of particle injection 
into the acceleration process and of energetic particle escape upstream and/or downstream 
of the shock. 

In this method, particle transport is described as a stochastic process. Particles 
move in small time steps, as their local plasma frame momenta are 'scattered' at each step 
in a random walk process on a sphere in momentum space. The properties of the random 
walk are determined by the assumption of a certain particle mean free path (or diffusion 
coefficient), that statistically describes the interactions of particles with the stochastic mag- 
netic fields. By assuming that such a description is possible, Monte Carlo methods gets a 
speed advantage over the PIC simulations at the cost of relying on theoretical models of the 
diffusive properties of the plasma. On the one hand, these models can be rather advanced 
and successful, therefore making this approximation justified. On the other hand, ignoring 
the spatial structure of electromagnetic fields and replacing it with a statistical description 
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is the biggest simplification of tliis model. 

Acceleration of particles takes place naturally in this model, as long as a shocked 
flow is described. Some shock heated thermal particles are injected into the acceleration 
process when their history of random scatterings in the downstream region takes them back 
upstream. These particles gain energy and some continue to be accelerated in the first- 
order Fermi mechanism. This form of injection is generally called 'thermal leakage' and 
was first used in the context of DS A in [48j (see also [l2] ) . The number of particles that do 
this back-crossing, and the energy they gain, are determined only by the random particle 
histories; no parameterization of the injection process is made other than the assumption 
of the diffusion coefficient value at various particle energies. 

The nonlinearity of the problem is dealt with by employing an iterative scheme 
that ensures the conservation of mass, momentum, and energy fiuxes, thus producing a self- 
consistent solution for a steady-state, plane shock, with particle injection and acceleration 
coupled to the bulk plasma fiow modification. 

An important advantage of the Monte Carlo model is that it was shown to agree 
well with spacecraft observations of the Earth's bow shock [50\ 151] . interplanetary shocks 
[8], and with 1-D hybrid PIC simulations 07]. 

2.4 Objectives of this dissertation 

Hopefully, I have convinced the reader of the far-reaching impact of processes in 
shocks on many astrophysical objects. Considering the observations of supernova remnants 
that indicate the possibility of strong magnetic field amplifications at shocks in-situ, I would 
like to theoretically investigate the physics of shock acceleration in the presence of strong 
MHD turbulence generation by the accelerated particles. 

I favor the Monte Carlo approach to this problem, because of the growing complex- 
ity of the models of nonlinear DSA (which makes the analytic approach less productive), 
and because I would like to probe the aspects of NL-DSA that neither the analytic models, 
nor the PIC simulations have yet constrained. 

The questions that interest me include: 

• How efficient can magnetic field amplification be, considering the nonlinear effects in 
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the system? What impact does the generated strong magnetic turbulence have on the 
efficiency of particle acceleration? 

• How do these results depend on the model of turbulence generation and on the model 
for statistical description of particle transport? 

• What are the consequences of efficient MFA on the maximum momentum of acceler- 
ated particles and how do they impact the shock structure? 

• Highest energy particles must escape upstream of the shock. What are the properties 
of the escaping particles? 

• What does the shock precursor look like (scale, structure, processes)? 

• What is the qualitative and quantitative dependence of some observable parameters 
(i.e., effective magnetic field strength, shocked gas temperature, flow compression 
ratio, etc.) on the properties of the shock (shock speed, plasma density and magneti- 
zation, etc.)? 

Over the past 3 years that the work on this project was being done, we (I, under 
the guidance of Prof. Ellison, and with the help of Prof. Bykov's advice) have successfully 
developed the model and obtained results shedding light on most of these questions. Our 
results have appeared in several peer-reviewed journal publications. 

The purpose of this dissertation is to make a record of the process of developing and 
testing the Monte Carlo simulation of NL-DSA with MFA, and to exhibit and summarize 
our results that have been or will soon be presented in conferences and in the press. 
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Chapter 3 

Model 

In the present research, I used the Monte Carlo method developed by Ellison and 
co-workers to build a self-consistent model of shock acceleration of charged particles, now 
with efficient magnetic field amplification. I wrote the computer code realizing the Monte 
Carlo model from scratch, but making a full use of the formerly developed procedures. I also 
contributed some essential improvements to the Monte Carlo method, that were necessary 
for the implementation of magnetic turbulence amplification models. 

In this Chapter, I will discuss the model. In Section [3.11 1 will present the funda- 
mentals of the Monte Carlo simulation and the tests performed to confirm that my numerical 
model reproduces the known analytic results and conforms with the fundamental laws of 
physics. Section [22] will be devoted to the state-of-the-art models of magnetic turbulence 
amplification discussed these days in the astrophysical literature, and to the implementa- 
tion of these models in the Monte Carlo code. This part of the model is the essence of my 
research project. Another original contribution I made to the model in this project is the 
adaptation and incorporation of advanced particle transport techniques into the simulation, 
as discussed in Section 13.31 Finally, in Section 13.41 1 will discss the realization of parallel 
computing in the simulation. 
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3.1 Core Monte Carlo 

This section discusses the techniques that were used in the Monte Carlo simulation 
before the incorporation of the magnetic turbulence amplification. Most of them had been 
developed before the author of this dissertation began contributing to the model. However, 
the tests of the model demonstrated in this section were performed by the computer code 
written by me. 

3.1.1 Overview 

The simulation of nonlinear particle acceleration with the Monte Carlo particle 
transport starts by assuming an unmodified shocked flow [u{x < 0) = uq, u[x > 0) = U2\. 
One must also assume some scattering properties of the medium, i.e., assign a mean free 
path A(x,p) to the whole particle energy range, at every point in space. 

Then thermal particles are introduced far upstream, and the code propagates these 
particles until they cross the subshock at x = 0. Particle propagation is diffusive, according 
to the chosen mean free path A(x,p), and it is performed as described in Section 13.1.21 
Some of these particles will be advected downstream with the flow U2, and once the code 
finds any particle many diffusion lengths downstream of the shock, its propagation may 
be terminated. However, the downstream flow speed must by definition be smaller than 
the downstream speed of sound (i.e., a thermal particle speed), so a small fraction of the 
downstream thermal particles may, in the random walk process, find themselves upstream. 
If this happens to a particle, it is said to have been injected into the acceleration process and 
becomes a CR particle (as opposed to having been a thermal one). An injected particle is 
much more likely than a thermal one to cross the shock again and again, and eventually gain 
a relativistic energy in this procesil]. The action of the advection with the flow combined 
with particle diffusion in a non-uniform flow is described in Section 13.1.31 

The particles that have been injected will, due to the flow speed difference across 

the shock, find themselves moving at a high speed with respect to the plasma, at least at 

the speed v = uq — U2- This completes the first cycle of the Fermi-I process. In a short time, 

the accelerated particles will again be advected downstream, but having a greater energy, 

^Note that in the colhsionless shocks discussed in this dissertation, particles with superthermal energies 
do not lose energy in particle-particle collisions and can therefore be easily accelerated once they are injected. 
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they will find it easier to return upstream again and get accelerated a second time. As this 
process goes on, a few particles may achieve, in principle, unlimited energy. 

The acceleration process in reality must have an upper limit for particle energy. 
The two possible causes of such a limit are a) time- limited acceleration: the particles only 
gain as much energy as they can in the amount of time that the shock has been in existence, 
and b) size-limited acceleration: particles have enough time to achieve such high energies 
that their scattering length becomes comparable to the size of the accelerating system, 
and they escape. The Monte Carlo method can, in principle, model both situations by 
terminating the acceleration of every given particle either after it had been in the process 
for a certain amount of time, or once it has reached a certain boundary. I favor the second 
scenario, due to its inherent consistency with the assumption of a steady state solution, and 
throughout this work it will be assumed that the acceleration is size-limited. For the model 
it means assuming a boundary, let us call it the free escape boundary (FEB) located at 
a^FEB < 0, such that any particle that crosses this boundary while moving against the flow 
leaves the system forever. The mean free path of the accelerated particles generally increases 
with energy, therefore only the highest energy particles can reach the distant boundary at 
XPEB- This way, the acceleration has an upper energy determined by the location of the 
free escape boundary. In reality, the distance xfeb must be comparable to, or be a fraction 
of, the radius of the SNR shell shock. 

After all the accelerated particles either had been advected downstream, or had 
escaped through the free escape boundary, the Monte Carlo transport process finishes. 
During the transport, the model was calculating the contributions of the particles, both 
thermal and CRs, to the particle distribution at every point. This process is described in 
Section [3.1.41 The fluxes of mass, momentum, and energy of the particles (and the accom- 
panying magnetic fields) are the important moments of the particle distribution function, 
that can be calculated directly from particle trajectories. 'Balancing the books' after the 
Monte Carlo iteration, one may either find that the above mentioned fiuxes were constant 
throughout the shock, or, if the particle injection was efficient, one may find that the calcu- 
lated fiuxes deviated from their upstream value. The latter case means that the accelerated 
particles gained too much energy from the bulk flow to be just a small admixture. Indeed, 
the Fermi-I process powers particle acceleration by giving a fraction of the bulk flow energy 
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to the particle scattering in it. There is only a finite amount of energy available for the 
accelerated particles, and as soon as they borrow a significant amount of it, the energy- 
bearing flow must change. This is where the nonlinearity of efficient particle acceleration 
comes into play. 

In order to obtain a steady state solution consistent with the fundamental con- 
servation laws (i.e., conservation of mass, momentum and energy), the simulation invokes 
an iterative procedure. Using the calculated, non-equilibrium, fluxes of mass, momentum 
and energy, the simulation adjusts the flow speed in the precursor, u{x < 0), making it 
decrease slightly towards the subshock, as outlined in Section TS.l.Sl It will be referred to as 
precursor smoothing. The value u{x = —0) = ui is the flow speed just before the subshock. 
Then a second Monte Carlo iteration of particle propagation may be run. Thermal parti- 
cles are introduced far upstream, where the flow is yet unmodified, u{x) = uo, and allowed 
to propagate, get injected and accelerated. However, this time, the flow speed difference 
between the downstream region and the upstream region is smaller due to the slowing down 
of the incoming flow in the precursor. This means that particle acceleration will borrow 
less energy from the flow than in the previous iteration. 

Continuing the iterative process of particle propagation followed by the flow speed 
adjustment, one may obtain a self-consistent solution, in which particles gain just the right 
amount of energy from the flow, to conserve mass, momentum and energy of the sum of the 
slowed down upstream flow and the accelerated particle distribution. It turns out that, in 
the presence of highly relativistic particles that change the compressibility of the plasma, 
and due to the assumption of the upstream particle escape at xpEB, one must also adjust the 
downstream flow speed, U2, in order to conserve the fundamental fluxes across the subshock 
as well. 

The procedure outlined above is the method of solving the problem of nonlinear 
particle acceleration developed by Ellison and co-workers. 

In order to implement efficient magnetic field amplification and self-consistent 
particle transport, the author made adjustments to the procedure. The underlying theory 
and practical details are explained in the following two Sections, 13.21 and 13. 3[ 

First of all, in addition to the flow speed, u{x) being an unknown function derived 
by the iterative procedure, we now seek for turbulence spectrum, W{x,k), in a similar 
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way (see Section l3.2p . Starting with some initial guess for W{x,k), the simulation runs 
the diffusion module, that analyzes the spectrum W{x, k) and calculates the corresponding 
mean free path at all energies, X{x,p) (Section 13.31 is devoted to this calculation). Then 
the Monte Carlo transport module is executed, that simulates particle acceleration in the 
given structure of the shock with the scattering properties determined by A. After that, 
collecting the information about mass, momentum, and energy fluxes, the model estimates 
the smoothing of the precursor required for the next iteration. Additionally, it collects 
the information about particle streaming during the acceleration process [i.e., the diffusive 
current of CRs, jd{x) or the CR pressure, Pcr(a;,p)]. Using this information, the code runs 
the magnetic field amplification module that calculates the turbulence generation by the 
particle streaming and adjusts the turbulence spectrum for the next iteration, W{x,k). 
The latter is then used to improve the guess on the particle mean free path, X{x,p), and 
the next Monte Carlo transport iteration starts. Similarly, this process continues until 
a self-consistent solution is derived, one that preserves the mass, momentum and energy 
fluxes, and in which particle acceleration produces the spectrum of accelerated particles 
that generates precisely the magnetic turbulence spectrum used for simulating the particle 
acceleration. 

3.1.2 Particle propagation, pitch angle scattering 
Theory 

Propagation of particles is performed using the methods developed and presented 
in [l9]. The bottom line of the reasoning provided in this work is the following procedure. 
If the particle has a mean free path A and a corresponding collision time tc = A/v, where v 
is the particle speed, then the scheme allows the particle to travel a finite time. At <C tc, 
in a straight line, and then rotates the particle's momentum. The momentum is rotated by 
an angle 69, which is chosen randomly as follows: 



cos 69 = I — {1 — cos A^max) • 



(3.1) 



Here X is a random number with a uniform distribution between and 1, and 




(3.2) 
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is the maximal scattering angle. The scheme works consistently when this angle is small. 
After choosing the polar angle of scattering, 56, one must choose the azimuthal direction of 
scattering, S(j). Assuming the scattering is isotropic, 

6(1) = 2TTy - TT, (3.3) 

where y is another random number uniformly distributed between and 1. To rotate the 
momentum, let's define spherical coordinates in the momentum space so that the azimuthal 
angle, 6, is measured from the positive px axis. Given the spherical angle of the original 
momentum, 9, such that cos^ = Px/p, we can calculate the spherical angle of the scattered 
momentum, 6', from: 

cos 6' = cos 9 cos 59 + sin 9 sin 59 cos Sep (3.4) 
The spherical angles (j) a-nd (p' don't matter in our 1-dimensional, axially-symmetric model. 




Figure 3.1: Pitch angle scattering diagram. 



The diagram illustrating the process is shown in Figure [3Tl The initial momentum 
vector p and the final (scattered) momentum vector p' are shown with thick arrows, and 
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the dotted lines crossing at the tails of p and p' show the maximal scattering cone with the 
half-opening angle of A^max- The vector p' is obtained from p by rotating the latter by a 
random angle 50 (such that {) < 69 < A^max) about the tail of p, and then by turning it by 
a random angle 5(j) (such that < (50 < 2ti) about the axis of p. 

Note that the term 'pitch angle' in plasma physics means the angle between the 
magnetic field vector and the particle momentum vector. This term was used in the name of 
this procedure because initially it was assumed that a uniform magnetic field, Bq, dominates 
the magnetization, with small fluctuations of this field providing the scattering. In the case 
of strong turbulence, however, the uniform pre-existing field will be overwhelmed by the self- 
generated fluctuations /S.B ^ i?o, so the term 'pitch angle scattering' is actually a misnomer. 
By the pitch angle in this scheme I simply mean the angle between the momentum vector 
p and the positive direction of the x-axis (the direction of upstream plasma flow) . 

Tests 

We illustrate the trajectories of the particles subject to the pitch angle scattering 
in Figure 13.21 For about 10 particles injected at some position at t = 0, the program 
recorded their deviations from the original position. Ax, versus time, t. It can be seen 
that the particles frequently change the direction of motion, which results in a stochastic 
transport. In the plot, the time t is normalized to the collision time, tcj which is defined as 
the ratio of the mean free path. A, to the particle speed, v. 

In order to verify that the properties of this particle transport correspond to dif- 
fusion, I plotted in Figure [3?3] the mean square of deviation of particle coordinates, (Ax^), 
as a function of time. The random walk process (i.e., such particle propagation that every 
step has the length A and is taken in a random direction) is a well-known textbook problem, 
and the solution predicts that the displacement \/ (Ar^) = \/iVA, where N is the number of 
the steps, and coordinate = x^ + + z'^. Therefore (|Ax^) = Dt, where D = v\/3, and 
t = Ntc- This dependence is shown in Figure [33] with the dashed line labelled "Theory", 
while the result of the Monte Carlo pitch-angle scattering simulation, averaged over 1000 
particles, is shown with the solid line. 

One can clearly see that there is a linear dependence of the mean square displace- 
ment, ^Ax^), on time, t. However, the slope of the solid line is some 20% greater than in 
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Figure 3.2: Particle trajectories calculated in the Monte Carlo code. 
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the random walk theory. This discrepancy is due to the definition of the mean free path 
made in [l9] in the derivation of the pitch angle scattering scheme. This definition is slightly 
different from the definition of the MFP in the random-walk model, where particles make 
discrete steps of length A in a completely random direction. However, a 20% difference in 
the value of the diffusion coefficient is a small factor compared to the uncertainties in the 
models of particle transport, and one can conclude that both models (random walk and 
pitch-angle scattering) describe the process of particle diffusion reasonably well. 

It must be explained why the authors of the Monte Carlo transport simulation 
chose the pitch angle scattering scheme instead of a simpler version of a random walk process 
with the step size being equal to the particle mean free path. The solution may, in principle, 
have sharp gradients of fiow speed, magnetic turbulence, or particle pressure, in which case 
the momentum distribution of some particles will be anisotropic. For example, the thermal 
particles first crossing the subshock, and finding themselves downstream, initially have a 
strongly anisotropic distribution of momentum. The random walk process assumes isotropy 
of particles in p-space, which may limit the physical veracity of the model, while pitch angle 
scattering accounts for the anisotropics exactly. In this sense, the Monte Carlo model 
goes one step beyond the diffusion approximation adopted in the simpler analytic models. 
This allows us to calculate the particle injection rate in the thermal leakage model, where 
the downstream thermal particles get injected by returning upstream in the process of 
their stochastic propagation, while some simpler analytic models resort to introducing an 
additional free parameter to fix particle injection rate (e.g., [22]). 

3.1.3 Motion of the scattering medium 
Theory 

With the diffusive transport of particles developed and tested, one needs to incor- 
porate the advection of the plasma. The paradigm of the Monte Carlo model (as well as of 
other approximations of the NL DSA problem) is that there is a bulk fiow of thermal plasma 
with respect to the shock, with magnetic fields frozen into it. Therefore, the scattering is 
elastic and isotropic in the plasma reference frame (elastic, because the scattering represents 
the action of the magnetic force, which doesn't perform work on the particle, and isotropic 
due to the assumption of strong developed turbulence adopted in the model; in some cases. 
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anisotropic scattering is a better approximation) . If there is a spatial variation in the speed 
of the bulk flow, then two subsequent particle scatterings may take place in two different 
reference frames, which may increase the energy of the particle (resulting eventually in the 
Fermi-I process). 

It is important to properly account for all effects of special relativity in order to 
model relativistic particle propagation. In order to model particle propagation in the moving 
plasma along with the diffusive transport, let us assume the following process. Suppose the 
particle is at the location x with the local plasma flow speed u{x). The code will allow this 
particle to travel a certain time in the plasma frame, and accordingly change its coordinate 
in the shock frame. Therefore, one may introduce an effective speed 

+ U{x) 

Vx, eflf = ; ■• (3.5) 

/ u^jx) 

Here Vx is the x-component of the particle's velocity in the plasma frame, and the physical 
meaning of eS is that, given a time At in the plasma frame, the particle's displacement 
in the shock frame will be Ax = Vx^ eflf^^- Then the program chooses the time the particle 
will be allowed to travel in the plasma frame. At, so that the resulting A^max will be small 
enough, as given by equation (j3.2p . After that, we change the particle's coordinate in the 
shock frame by Ax = Vx, eff At and proceed with the pitch angle scattering routine. 

An important new aspect of this process is that in a model with efficient magnetic 
field amplification, one expects large gradients of the mean free path as well as of u{x). 
This means that when Ax is large enough for the particle to cross one or more numerical 
grid planes at which all the physical quantities are defined (see Section I3.1.4p , care must 
be taken to account for the changing properties of the medium. The simulation deals with 
this in the following manner. I choose a single numerical value of the angle AQmax for all 
particles throughout the simulation. Then the code allows each particle to propagate just 
enough time, so that the 'accumulated' maximal scattering angle is exactly AOmax- If a 
particle is to cross several grid planes in the course of this time, then I define a cumulative 
maximal scattering angle as 

AeLx,cmi = 65]^, (3.6) 

^ TcoU, i 

where the summation index i runs over all the spatial bins that the particle had crossed, 
each bin with a different value of TcoH, and Atj is the amount of time that this particle 
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spent in bin i. Only after AOmax, cmi = ^©max! does the Monte Carlo routine scatter the 
particle. This makes the results independent of the choice of grid plane locations, as long 
as the separation between them is small enough. 

Test 

Let us perform a test of the advection superimposed on diffusion in the Monte 
Carlo code. With a bulk plasma flow of speed uq set up, I introduce the test particles at 
X = with their plasma frame velocity in the positive x direction, and let them propagate. 
Particles are assumed to move according to the pitch angle scattering scheme described 
in the previous section; between the scatterings, the motion of the particles is ballistic in 
the plasma frame, which moves at the speed uq with respect to the stationary frame. The 
scatterings are isotropic and elastic in the plasma frame, with the corresponding momenta 
in the stationary frame Lorentz transformed. 
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Figure 3.4: Advection with diffusion. 

In Figure 13.41 I illustrate the results of the test. I introduced 3 particles with 
different speeds in the plasma frame: a slow particle with v = uq/5, a moderately fast 
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particle with v = uq, and a fast particle with v = 5uo (here and below, the letter u will 
denote the speed of the flow, and the letter v will refer to the speed of a particle, measured 
either in the stationary, or in the plasma reference frame). The solid line represents the 
motion of the first, the dashed line - of the second, and the dash-dotted line - of the third 
particle, respectively. The spatial coordinate, x, is measured in the units of rgo, the latter 
being the mean free path of the particle with speed v = uq, and time is measured in units 
of tgo = rgo/uQ. It was assumed that the mean free path is proportional to the speed of the 
particle in the plasma frame. One can see from the Figure [33] that the slow particle moves 
almost synchronously with the flow (the solid line is very close to x = uot). The faster 
particle, with v = uq, on the average moves along with the plasma, but sometimes slower, 
and sometimes faster, because its x-component of velocity in the stationary frame, Vx -\-uq, 
varies between and 2uq, depending on the orientation of the particle's momentum. The 
fastest particle, v > uq, can move backwards, as the dash-dotted line shows. However, its 
motion is still affected by the flow, shifting the average location at the advection speed uq. 

Because the purpose of the model under development is to model nonlinear shock 
acceleration, where the flow speed u{x) can vary with distance, smoothly (in the shock 
precursor) or discontinuously (across the subshock), one must see how the code treats such 
varying flow speeds. Let us study three cases: a uniform flow speed, a smoothly varying 
flow speed, and a discontinuity in the flow speed (representing a shock). In a separate test, I 
will introduce a slow particle {v = O.Suq) at x = —100 rgo, where u{x) = uq, and trace it as 
it propagates in the flow. The code will record the positions of the particle, measured in rgo 
(the latter is, again, the mean free path of a particle with v = uq, and the mean free path 
for any other particle energy is proportional to the particle speed), and the x-components 
of the particle's velocity in the plasma frame, Vx- 

Figure [331 shows the results for these three cases. In the first case (top panel), 
with constant flow speed, Vx varies with time in the range — O.Sno < Vx < +0.3tio, with 
the average (vx) = 0, but the dispersion (^v'^'j remains constant. This case is similar to the 
situation studied in the previous test. In the second case (middle panel of Figure [33|) . where 
the flow speed linearly drops from u{x = —100 rgo) = uq to u{x = 0) = uq/10 and then 
remains constant, we see a different behavior. The average particle motion is still locked 
with the plasma {{vx) = 0), but the dispersion (v^) increases as n(x) drops. This means 



41 




42 



that the particle's energy in the plasma frame grows. Such energizing of the particles is, in 
fact, the adiabatic heating of a gas put in a slowly shrinking volume (see Appendix B of 
[123| or Section 1.2 in [126J). In the third case (bottom panel of Figure [33|) . where the flow 
speed is constant at x < 0, but drops abruptly from u{x = —0) = uq to u{x = +0) = no/10, 
similarly to what it looks like in a shock, the particle is energized significantly at the shock 
crossing. I, actually, had to introduce a reflecting boundary at x = that doesn't allow 
the particles to cross the shock backwards, from x > to x < 0, in order to show a concise 
plot. Such crossing back becomes possible because the speed of the particle in the plasma 
frame is greater than the shock speed. 

3.1.4 Calculating particle distribution and its moments 

To calculate the particle distribution in the simulation, the simulation registers 
particles crossing certain locations, that we hereafter refer to as 'grid planes', because these 
locations also define the spatial grid, at the nodes of which all the quantities: u{x), W{x, k), 
etc., are defined. One may think of this as detection of particles by imaginary detectors 
placed at discrete locations upstream and downstream of the shock. This calculation may 
seem a bit tricky because each crossing of a detector by a particle contributes not to the 
density, but to the flux of particles, so in order to extract the particle distribution infor- 
mation, one needs to properly weight the detected information. However, this weighting 
is a standard procedure for simulations like ours, and below I demonstrate the reasoning 
leading to it and examples of the scheme at work. 

The Monte Carlo simulation does not populate the whole space with particles; 
instead, it introduces Np particles upstream and propagates them one by one, until each 
leaves the system; yet, it is simulating a steady state solution with this process. This means 
that if one wants to calculate the particle distribution function /(x,p) and its moments 
(momentum and energy fluxes) at some spatial locations, one must collect the information 
about the particles in such a way that all the data collected in the course of one iteration 
(i.e., during the propagation of all the Np simulation particles) represents the information 
that would be collected by particle detectors placed in the plasma in a unit time (e.g., in 
one second). 

The simulation only registers the particles' contribution to f{x, p) when they cross 
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one of the 'detectors'. Consider an infinite plane detector in a spatially uniform plasma with 
no bulk motion {uq = 0) and an isotropic distribution of particle momenta; assume also that 
all particles have the same speed v. Clearly, in a unit time such a detector will register more 
particles incident normally onto it than tangentially to its plane. It is easy to understand 
that the number of detected particles with a certain Vx is proportional to P{vx) oc \vx\, 
where Vx is the x-component of the particle's velocity in the rest frame of the detector and 
P{vx) is the probability density of its detection in a unit time. 

When a 'detector' in the simulation registers a particle with small \vx\, we must 
interpret it as that there are many similar particles at this location, but, because of an 
unfavorable direction of motion, only a few reach this parallel-plane detector in a unit time 
that one iteration represents. 

Quantitatively, if one wants to calculate the number density of particles at the 
location of the detector, the code must compute the following sum: 

3 3 

where i is the number of the grid plane, and Xj - its coordinate, n{xi) is the number 
density of particles at Xj, and the summation index j runs over all the events of a particle 
crossing this detector in the course of an iteration. The weight Wj is the statistical weight 
of the j-th event, and Wp is the statistical weight of the particle participating in the event. 
The statistical weight of a particle, if Np particles are introduced upstream, representing a 
plasma density ng, is simply Wp = rio/Np. The statistical weight of the event is expressed by 
the ratio \uo/vx,j\, where the denominator Vx^j is the x-componcnt of the particle velocity 
measured in the rest frame of the detector, at the j-th event, and uq acts as the appropriate 
normalization factor. 

If one wants to calculate the particle distribution function, i.e., the number of 
particles in a unit phase space volume dxdydzdpxdpydpz, then the summation must be 
restricted to the events corresponding to that phase space volume. In practice, one may be 
interested in the distribution function /(x,p) such that 

J f{xi, p) d^p = J f{xi, p)p^dpdn = n{x), (3.8) 

where dfl represents the infinitesimal spherical angle corresponding to the momentum space 
volume d^p = dpxdpydpz- In the simulation, given a phase space binned so that Apk is the 



(3.7) 



44 



width of the k-th momentum bin centered at the momentum value pk, and averaging over 
the angles, one gets 



f{xi,pk) = ^ I f{xi,p) dn 
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where dVL is the differential of the solid angle in p-space, and the index j runs over all 
particles whose momentum falls into the k-ih momentum bin. Indeed, with f{x,p) defined 
this way, 

f{xi,p)d^p = f{xi,p)p'^dpdn = 
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as expected (in the summation, the index k runs over all the momentum bins defined in the 
simulation, and in the last summation, j runs over all the events of particle crossing of the 
detector). In the future, the plots showing f{x,p) actually show f(x,p), but I omit the bar 
representing the angular averaging for simplicity. 

Calculating the angle-averaged distribution function may be informative, but for 
practical purposes one needs the moments of the distribution function (i.e., the mass, mo- 
mentum and density fluxes) that require the information about the angular dependence. 
From the above reasoning (namely. Equation I3.10p one may conclude that the correspon- 
dence between the integration over the momentum space and summation over particle 
detection events is as follows: 
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f{x,p)d^p= 

PjGApk 

Therefore, for any function of coordinate and momenta Ai{x,p), its expectation value may 
be calculated directly in the simulation as 
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J M{xi,p)f{xi,p) d^p = YM{x,pj] 
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Namely, for Ai{x,p) = 1 one finds the particle density, for Ai{x,p) = mpVx - the mass 
flux in the a;-direction, for Ai{x, p) = PxVx — the flux of the x-component of momentum in 
the x-direction, and for A^(a;,p) = K{p)vx - the energy flux in the x-direction {K{p) is the 
relativistic kinetic energy corresponding to the momentum p = |p|). The expectation value 
of the function M{x,p) = ev is the diffusive current jd{x), and M.{x,p) = ^pv in the case 
of isotropic momentum distribution gives the pressure. 
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(3.13) 
(3.14) 
(3.15) 
(3.16) 
(3.17) 
(3.18) 
(3.19) 



etc. In the last two equations, the integration was limited to the thermal or to CR particles 
only. In this way I define the thermal pressure and the CR pressure. I must remind the 
reader here that a peculiarity of the approach I adopted is that in order to separate the CR 
particles from the thermal ones, the code uses their history, and not their energy. By my 
definition, a thermal particle is one that had been introduced into the simulation upstream 
with a random thermal energy and that may have crossed the subshock going downstream, 
but has never crossed it back. Once a particle crosses the subshock (the coordinate a; = 0, 
to be more precise) in the upstream direction, it by my definition is injected and becomes 
a CR particle. 

Let us note here that, assuming that the particle distribution is isotropic, the 
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fluxes of energy and momentum can be expressed via gas pressure as 

j p^v^f{x,p)d^p = p{x)u^{x) + Pp{x), (3.20) 

Kv,J{x,p)d^p = ^p{x)u^{x) +wp{x)u{x) , (3.21) 

where Pp{x) is the pressure and Wp{x) is the enthalpy of the particles. These are well known 
results of the kinetic theory of gases. 

3.1.5 Introducing particles into the simulation 
Theory 

In order to start the Monte Carlo simulation of particle transport, we must intro- 
duce thermal particles far upstream, or close to the subshock (the latter method is described 
in Appendix B of [123 J). This task has two components: generating a particle population 
with the proper energy distribution, and choosing the appropriate angular distribution for 
these particles. 

We normally assume that the unshocked plasma is thermal and has a certain 
temperature Tq (a typical cold interstellar plasma has Tq ~ 10"^ K). In order to generate a 
thermal population, we randomly choose the momentum of every particle introduced into 
the simulation so that the resulting distribution is Maxwellian: 

/(p) = no( ^— ^ ^ exp ( . (3.22) 

\2TTmkBTo J \ ImKBTo J 

In order to accomplish this, consider the function F{p) such that F{p)Ap is the fraction of 
the thermal particles with momenta in the interval [p — Ap/2;p + Ap/2], i.e., 

4V/(p) 4 f 1 y/^ 2 / p^ \ 

F[P) = = ^ ^ , ^ p exp - . (3.23) 

no ^/^T \2mkBToJ \ 2mkBTo J 

By substitution 



using the identity G{y)dy = F{p)dp, we find that y has the following distribution function: 

G{y) = -^yh-y (3.25) 
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(here G{y)Ay is the fraction of particles with momenta corresponding to the interval [y — 
Ay/2; y + Ay/2]). The latter is a Gamma distribution with parameter a = 3/2. To generate 
a quantity y with the above distribution, one can use the following recipe: 

y = ^Z^-lnA:. (3.26) 

Here X is a random deviate with a uniform distribution in (0; 1], and Z is a deviate the 
the normal (Gaussian) distribution with the mean equal to and the dispersion equal to 1. 
The above method was adopted from |37j . 

Test 

I tested the implementation of this procedure, and verified the distribution function 
in the Monte Carlo simulation, and the test results are shown in Figure \3M 
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Figure 3.6: Generation and detection of thermal particle population. 



The code introduced Np = 10^ particles into the simulation, assuming that their 



temperature is Tq = 10^ K, and the density no = 0.3 cm ^. The thick line shows the 
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desired Maxwellian distribution calculated according to Equation (j3.22|) . The thin line 
shows the result of the detection of the introduced particles at some grid plane, as described 
in Section [3. 1.41 according to Equation (j3.9p . The local deviations of the Monte Carlo result 
from the theoretical curve are statistical fluctuations, and these decrease for greater number 
of particles. Otherwise, the match is excellent, demonstrating the correct implementation 
of the introduction and the detection of particles in my code. 



of concern for a simulation like ours, because it determines the rate of particle injection 
into the acceleration process. When the simulation introduces particles at the coordinate 
X, it is replacing the dynamics of these particles upstream of x with an analytic description, 
consequently it must distribute particles in p-space at x the way they would be distributed 
having traveled from far upstream and reaching x for the first time. This is equivalent to 
calculating a p-space distribution of particles incident on a fully absorbing boundary at x 
after scattering in a non-uniform flow u{x). This is easy to do analytically if all particles 
have a plasma frame speed v less than the flow speed u{x) (because then all particles 
crossing position x do it for the first and the last time), and fairly complicated otherwise 
(see Appendix B). Let us assume v < u{x) in further reasoning, which is justified as long 
as the local sonic Mach number at the introduction position is large. 



introduced particle's momentum, p. But how do we choose its direction, identified by 
the angle /x such that cos /x = Px/p'^ As was stated earlier, it is assumed that the angular 
distribution of momenta of the introduced thermal particles is isotropic in the plasma frame, 
and there is an overall drift speed u superimposed over this motion in the plasma. Therefore, 
it may seem natural (but is incorrect) to just choose px isotropically in the plasma frame, 
and then transform them into the shock frame. The correct solution must account for the 
fact that when these particles cross a grid plane, their flux must be 'fiux-weighted' as seen in 
Equation (|3.7|) . because the number of particles arriving at x in a unit time is proportional 
to the cosine of the angle that their shock frame velocity Vgf makes with the x-axis. This 
can be done by assuming a probability density of of Vsf^ x as 



The angular distribution of momenta of the introduced particles is a major issue 



The problem is now reduced to the following. We know how to designate an 




(3.27) 
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Here fmin = u — v, Umax = u + V, and V is the particle speed in the plasma frame chosen 
using a random number generator according to (j3.26p . The constant A can be found from 
the normalization condition 



as 



A = —, (3.29) 

V — V ■ 



SO 

2Usf , X 



i^(t^sf,x) = < 



'^min ^ ^sf , X ^ ^11 



71^ — V ■ 



(3.30) 



0, otherwise 
In order to generate such a particle distribution, we use a random number Z uniformly 
distributed between and 1 and calculate Wsf, x as a function of Z. The identity 

F«)d< = H{z')dz', (3.31) 

and substitution of the uniform distribution 

. ,x f 1' < ^' < 1> .X 

i/(/) = <^ (3.32) 
I 0, otherwise 

lead to: 



F{v'Jdv'^ = J dz' = z, (3.33) 



which is where we can derive Ugf, x(^) from. Substituting ()3.30p into (j3.33p . we get 



^sf, X = Y «ax - vl,iJZ + i;^;^. (3.34) 

Note that, if we did not account for the flux weighting and just prepared an 
isotropic distribution of particles in the plasma frame and then transformed it into the 
shock frame, then instead of the prescription (j3.34p we would use 

Vs{^^ = {2Z -l)v + u, (3.35) 

which is incorrect, as I show below. 

If properly implemented, the introduced particle population as detected by the 
grid plane detectors should have a uniform distribution of /x in the plasma frame, g{fi) = 
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1/2. This corresponds to isotropy of /(p). Let us perform two tests to confirm that the 
procedure (j3.34p gives the correct particle distribution isotropic in the plasma frame. 

In the first test, let us introduce Np = 10^ particles into the simulation with a 
supersonic flow (sonic Mach number = 2.5) using the incorrect recipe (|3.35|) . At the 
grid plane very close to the introduction position, the model will measure the angular dis- 
tribution of particles, g{^). Several (approximately 20) diffusion lengths downstream of the 
introduction position, it will measure g{^) again. By the time of the second measurement, 
the particles must have scattered enough to assume an isotropic velocity distribution in the 
plasma frame corresponding to g{^) = 1/2. The results of the test are shown in Figure [321 

As expected, the thick line in Figure [3771 is g{ii) = 1/2, meaning that after many 
scatterings, particles isotropize their velocities, and also confirming that I implemented 
correctly the calculation of particle distribution and the pitch angle scattering routine. At 
the same time, the tilt of the thin line tells us that the introduced particle distribution was 
not isotropic in the plasma frame, thus the recipe (j3.35p is incorrect. 

In the second test, Np = 10^ particles will be introduced, now using the recipe 
(j3.34p to choose their angular distribution. Let the model measure the angular distribution 
immediately, and some distance downstream of the introduction position. The results are 
shown in Figure 13.81 The two angular distributions match exactly, which means that the 
introduced angular distribution was isotropic, and remained such after many scatterings. 

These tests conclude the verification of particle propagation methods of Monte 
Carlo, and now we can get to testing the particle acceleration properties. 

3.1.6 Test particle case of DSA 

A crucial test that the nonlinear simulation of DSA must pass is the production 
of a power-law spectrum of accelerated particle in the test particle case. Similarly to the 
solution shown in Section ri.5[ a shock characterized by a sharp transition at x = 0, in which 
the flow speed jumps from the uq to U2, with a compression ratio r = uq/u2 > 1 and particle 
injection taking place at the shock via thermal leakage, must produce a power-law spectrum 
of accelerated particles with the power law index s given by Equation (|1.23p . Verifying the 
power law index is a strong argument for the correctness of the model. 

We ran 3 simulations with a discontinuous flow speed. In each simulation, uq = 
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Figure 3.8: Test of introduced particle distribution isotropy. 
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10^ km s~^, To = 10^ K, uq = 0.3 cm~^, and to define the mean free path, the model assumed 
Bohm scattering in a magnetic field Bq = 3 ■ 10~^ G (see Section I3.3.ip . The difference 
between the 3 models was the compression ratio, r. I chose r = 7.0 for the first model, 
r = 4.0 for the second and r = 3.0 for the third. Note that I did not require consistency 
of these compression ratios with the laws of hydrodynamics, and was only interested in 
particle acceleration and its properties. Namely, according to Equation (jl.23p . one expects 
to get power-law distributions of accelerated particles downstream with the indices s = 3.5, 
s = 4.0 and s = 4.5, for the first, second and the third model, respectively. In these runs 
I took advantage of the procedure of particle introduction developed in Section 13.1.51 and 
introduced particles close to the shock instead of far upstream. This allowed us to speed up 
the calculation significantly, and thus I put the free escape boundary rather far upstream, 

at XFEB = —10^ TgO. 




Figure 3.9: Test particle case of DSA. 

The results of the test - the particle distributions measured in the shock rest 
frame - are shown in Figure 13.91 The solid lines are the particle distributions in the shock 
frame measured in the downstream region. The j;-axis shows proton momentum in units of 
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nipC, and the y-axis - the distribution function f{p) multiphed by p'^ for convenience. This 
multipUcation factor makes the r = 4.0 case with the corresponding s = 4.0 appear as a 
horizontal hne, which is a desirable feature of this plot. In the future, all particle distribution 
functions shown will have this multiplication factor. The thin solid line shows the detected 
distribution function for r = 7.0, the medium thickness solid line - for r = 4.0 and the thick 
solid line - for r = 3.0 case. The dashed lines are added for a comparison. They have the 
slopes predicted by Equation (jl.23p . and for the correct result, the dashed lines must be 
parallel to the high-energy parts of the distribution functions, which is certainly the case in 
the presented results. 

I would also like to illustrate the physical process that leads to particle acceleration 
in these tests. In a separate simulation similar to the r = 7.0 case, but with a free escape 
boundary close to the shock, at xfeb = —50 rgo, the code traced several particles, and I 
show their motion in the phase space in Figure 13.101 

The top panel of Figure [3.101 shows the shocked flow speed. Upstream, at x < 0, 
the flow is uniform and fast, uq = 10'* km s~^, and at the shock located at a; = 0, the flow 
speed drops down to U2 = uo/r. This flow speed is measured in the system where the shock 
is stationary, so the flow of matter is directed from the left end to the right end of the plot, 
and the shock is directed to the left. 

The second panel shows a particle introduced far upstream, that crossed the shock, 
got heated, but never returned upstream and escaped trapped in the downstream flow (the 
particle's motion is from the left to the right end of the plot). This is what happens 
to all particles in collisional shocks that are usually observed on Earth, and no particle 
acceleration occurs. 

The third panel shows a particle initially introduced far upstream as a thermal 
particle, but its random motion in the stochastic magnetic fields, induced by the pitch angle 
scattering, lead the particle back upstream, and it crossed the shock from the x > region 
into the x < region. In several such subsequent crossings, the particle gained energy due to 
the flow speed difference across x = (the trajectory moved up). Eventually this particle's 
energy became so large that it was able to protrude quite far upstream, to x ~ —10 rgo, 
but eventually it was advected downstream with the flow. 

The fourth (bottom) panel in Figure 13.101 shows the phase space trajectory of a 
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Figure 3.10: Particle trajectories in DSA. 
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'lucky' particle that not only got injected (crossed the shock against the flow), but gained 
enough energy to find itself very far upstream, to the left of the free escape boundary 
located at a; = —50 rgo- There such particles were assumed to escape from the system 
in the upstream direction - this is how the finite size of the accelerator is modeled in the 
Monte Carlo simulation. The energy of this particle at the moment of escape was close to the 
maximum achievable particle energy in this shock. The particle's speed greatly exceeded the 
shock speed at this moment, with the particle actually being mildly relativistic p ^ 0.6 iripC. 
In a model representing a real SNR shock, the free escape boundary would have been 
located much farther upstream, and particles would be able to make many more crossings 
before escaping upstream, and the maximum particle momenta could be ultrarelativistic, 
p 10® nipC. 

3.1.7 Shock compression ratio in nonlineeir DSA 

Theory 

In a standard steady state hydrodynamic shock, the relationship between the 
pre-shock macroscopic quantities: uq, uq, Tq and their post-shock values: U2, n2, T2 is 
determined by the sonic Mach number of the shock, and can be derived from the Rankine- 
Hugoniot equations: 

P2U2 = poUu, (3.36) 
P2UI + P2 = Poul + Po, (3.37) 
-P2ul + W2U2 = -^PoUo + wquq, (3.38) 

which express the conservation of mass, momentum and energy fluxes, respectively. Here 
P is the gas pressure, and w - its enthalpy, w = e + P, and the internal energy e of the 
gas is proportional to the pressure P. For an adiabatic gas with the ratio of specific heats 
7, one can write w = "yP/i^ — 1), and the upstream gas pressure Pq can be related to 
the sonic Mach number Mg as Mg = (uq/cs)"^ = pqUq/^jPo). This leads to a solution of 
the Rankine-Hugoniot equations, relating the pre-shock and the post-shock flow speed and 
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temperature. That solution is called the Hugoniot adiabat: 

uq 7 + 1 



U2 7 + 2/M2 - 1 ' 

_ (27M2-(7-1))(2/M2 + (7-1)) 



(3.39) 

2 , (3.40) 



To (7 + 1 

where uq/u2 = P2I Po = ^tot is the compression raticj^. 

In a nonlinear shock with magnetic field amplification, the situation is complicated 
by the contributions of cosmic ray pressure and magnetic turbulence pressure, and by par- 
ticle escape far upstream. The procedure of the search of the self-consistent compression 
ratio, developed in [l3], is based on the requirement that in a steady-state system, mass, 
momentum and energy fluxes must be constant in space. I generalized this procedure for the 
problem of magnetic field amplification by including the contributions of magnetic turbu- 
lence to the conservation relations, and developed an iterative procedure for an automated 
search of the self-consistent rtot • Consider the conservation relations 

p(x)u{x) = Polio (3.41) 
$p(x) = $Po, (3.42) 
$i?(x) + Qesc(x) = $E0. (3.43) 

Here p and u are the mass density and the flow speed, $p(x) is the flux of the x-component 
of momentum in the x-direction including the contributions from particles and turbulence, 
and <I>po is the far upstream value of momentum flux, i.e., 

^PO = Powo + ^tho + ^wo • (3.44) 

The quantity <I>p is defined as 

^P{x) = j p,vJ{x,p)d^p + P^ix), (3.45) 

where Px and Vx are the x-components of momentum and velocity of particles, and /(x, p) 
is their distribution function, all measured in the shock frame. The quantity <I'_e(x) is the 
energy fiux of particles and turbulence in the x-direction, Qqsc 

is the energy fiux of escaping 



^Hereafter let us replace the notation of the total shock compression ratio. Instead of r, we will now 
denote it as not = uo/u2, to distinguish it from the subshock compression r-sub = uxju^. 
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particles at the FEbI^I and the far upstream value of the energy flux is 

^Eo = IpouI + Ptho^o + Fwo- (3.46) 
z 7 ~ J- 

The quantity ^e{x) is defined as 

^e{x) = J Kv,f{x,p)d^p + F^{x), (3.47) 

K being the kinetic energy of a particle with momentum p measured in the shock frame, 
and are the momentum and energy fluxes of the turbulence defined in Section 13.2.91 
Writing equations (I3.4ip . (13.421) and (I3.43P for a point downstream of the shock, 
sufficiently far from it that the distribution of particle momenta is isotropic, and the approx- 
imations (|3.20p and (|3.2ip are valid, and denoting the corresponding quantities by index 
'2', we get the equivalent of the Rankine-Hugoniot relations, that accounts for particle 
acceleration and escape, and for the presence of magnetic turbulence: 

P2U2 = PqUq, (3.48) 
P2ul + Pp2 + Pw2 = Poul + PpO + PwO = '^PO, (3.49) 

1 o Is 

-P2U2 + Wp2U2 + -Fw2 + Qesc = -^PQUq + ''^pO'U-O + -^wO = ^EQ. (3.50) 

The particle gas enthalpy Wp is Wp = tp+Pp, and the internal energy ep of gas is proportional 
to the pressure Pp. Introducing the quantity 7 so that ep = Pp/{^ — 1), one can write 

WpU = _ ^ ^ PpU (3.51) 

The value of 7 is averaged over the whole particle spectrum, and it ranges between 5/3 for 
a nonrelativistic and 4/3 for an ultra-relativistic gas. The local value of 7 can be easily 
calculated in our code from the particle distribution, along with Pp and €p, as 7 = l + Pp/cp. 
Similarly, one can define 5 = Fm/{uPw) and calculate a local value of 6 anywhere in the 
code in order to express 

F^ = 6- PyjU . (3.52) 

The value of 5 depends on the nature of the turbulence. For instance, in Alfvenic turbulence, 
one expects 5^3, [see Equation (j3.178p ]. 



^Particle escape at an upstream FEB also causes the mass and momentum fluxes to change but these 
changes are neghgible as long as uq ^ c (see [43]). 
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Substituting (j3.51|) and (|3.52|) into the above equations and introducing rtot = 
uq/u2., we can eliminate p2 using (j3.48p and Pp2 using (j3.49p . which ahows us to express 
from ()3.50p the quantity qcsc = Qcsc/^EQ as 

.esc = 1 + (3.53) 



where 



C 



A = (3.54) 



72-1' 

B = 1 + -^^ 2 2~' W-^^) 

72 - 1 V PoUq J PoUq 



270 Ppo ^ 2<5oP„ o 
70-1 po^o Po^^o 



C = (3.56) 



Note that po^o/-^pO = 7o^s > where 70 = 7 = 5/3 due to the absence of CRs far upstream. 
The pressure of stochastic magnetic fields PwO can be found from the spectrum of seed 
turbulence far upstream (see Section 13. 2p . 

The quantity gesc is readily available in the simulation after the end of any iteration. 
Comparing it to the value predicted by (j3.53p , one may evaluate the self-consistency of the 
solution and make the correction to rtot, if necessary, for further iterations. For making 
these corrections it is helpful to use in the simulation the inverse of ()3.53p . the physically 
relevant branch of which is 

2A 

= B-^B'^-4AC{l-q^.y ^'-''^ 
It is important to emphasize here that an iterative procedure is required to find 
the compression ratio rtot of a non-linear ly modified shock, because quantities qesc, Pw2 and 
72 depend on rtoti so (j3.57p only provides a practical way to perform the iterations. The 
code employs the following procedure: 

2A 

r(„ = (1 - ()rtot + C , (3.58) 

^ij_^i32-4^C(l-gesc) 

where (" is a small number (typically ( = 0.01 . . . 0.1). Here rtot is the compression ratio 
assumed for the last iteration, and r[^^ - the compression ratio chosen for the following 
iteration. The weighting using the parameter is chosen so that, when the deviation of 
the nonlinear structure of the shock from the self-consistent solution is large, this iterative 
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procedure would not overcompensate the discrepancy, which may lead to the breakdown 
of the model (for example, rtot < 1 is unphysical, and rtot too high may stall the particle 
transport procedure). 

Another complication that may arise with the procedure described by (|3.58|) is 
that instead of converging to a solution that satisfies the conservation relations (j3.48p . 
()3.49p and (13.501) . the procedure may find an attracting cycle around the self-consistent 
value of rtot- For example, rtot niay be too low in one iteration, underestimating particle 
acceleration efficiency, which leads to a shock with a high r[^^ predicted by (j3.58p . The latter, 
in turn, overestimates particle acceleration, leading to a low compensatory rtot prediction 
again. These cycles may be merely a mathematical consequence of our numerical model 
of particle accelerating shocks. On the other hand, it is conceivable that a real shock, 
instead of evolving into a steady-state system, may behave periodically or even chaotically, 
turning particle acceleration on and off. However, these effects are beyond the scope of the 
present research, because we look for the steady-state structure of collisionless shocks. The 
simulation avoids the attractors other than the self-consistent solution by randomizing the 
value of C, so that it varies between a finite value and in every iteration. This way, any 
attracting cycle that the system (j3.58p may have with a constant value of (" will eventually 
be broken, but if rtot is at its self-consistent value {r[^i ~ ftot)^ the randomization of C will 
not take the solution away from this point. The latter is expected as long as statistical 
fluctuations of quantities A, B, C and gesc keep r[^^ in the attracting domain. This only 
requires that a high enough number of particles is used in the Monte Carlo routine. 
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Test of the implementation 

We test the procedure of the iterative estimation of the compression ratio (j3.58p by 
confirming that it reproduces the solid predictions of the Hugoniot adiabat (j3.39p and (j3.40p . 
In 10 runs described below, the shocks propagate in a gas with density no = 0.3 cm~^, 
temperature Tq = 7.3 • 10^ K (corresponding to a sound speed Cs = 10 km s^^), and 
magnetic field Bq = 10~^ G determining the Bohm diffusion (this magnetic field is too 
small to influence the momentum and energy balance, therefore the Hugoniot adiabat should 
apply). I performed 9 runs, in which the flow speed, no, varied from 5 km s^^ to 1000 km s^^, 
corresponding to the sonic Mach number, Mg, varying from 0.5 to 100. In these simulations, 
I artificially eliminated particle acceleration by assuming that the subshock is fully reflective 
for particles trying to cross it from the downstream into the upstream region. Additionally, 
I ran simulation number 10, for which the subshock is assumed fully transparent, and 
particle acceleration occurs (limited by a free escape boundary at x = —80 Tgo). We start 
the simulations off by assuming a flow with the speed uq, which at the point x = abruptly 
slows down by 0.1% to U2 = uo/l-OOl. This tiny flow speed jump heats the particles a 
little, just enough to make the procedure (|3.58p start converging to a self-consistent rtot- 
Setting the parameter = 0.3 and randomizing it between and 0.3, in 30 iterations the 
simulation obtains the self-consistent value of rtot- I averaged the rtot prediction over the 
last 10 iterations out of 30, and showed the results in Table [3Tl Additionally, the simulation 
measured the downstream gas temperature, T2, by detecting thermal particle pressure Pth2 
and relating it to the temperature by the ideal gas law. 

Column 'Model' in Table 13.11 shows the number of the model, uq is the upstream 
flow speed, Mg ~ the corresponding Mach number, 'Accel.' shows whether the Fermi-I 



Table 3.1: Test of iterative search of the compression ratio, rtot 



Model 


uo, km/s 


Ms 


Accel. 


(not) HA 


(not) MO 




(T2/To)mo - MC 


1 


5 


0.5 


no 


1.0* 


1.022± 0.003 


1.0* 


0.99 ± 0.08 


2 


10 


1.0 


no 


1.0 


1.05 ± 0.01 


1.0 


1.02 ± 0.05 


3 


12 


1.2 


no 


1.30 


1.31 ± 0.02 


1.19 


1.19 ± 0.06 


4 


15 


1.5 


no 


1.71 


1.73 ± 0.01 


1.49 


1.48 ± 0.06 


5 


20 


2.0 


no 


2.29 


2.29 ± 0.01 


2.08 


2.0 ± 0.1 


6 


30 


3.0 


no 


3.00 


2.99 ± 0.01 


3.67 


3.7 ± 0.1 


7 


50 


5.0 


no 


3.57 


3.58 ± 0.02 


8.68 


8.4 ± 0.3 


8 


100 


10 


no 


3.88 


3.90 ± 0.01 


32.1 


31 ± 1 


9 


300 


30 


no 


3.99 


4.00 ± 0.01 


282 


290 ± 10 


10 


300 


30 


yes 


3.99** 


3.11 ± 0.09 


282** 


230 ± 10 
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acceleration was allowed, as described above. The columns '(rtot)HA' '(-^2/7o)ha' 
the predictions of the Hugoniot adiabat (|3.39p and p.40p corresponding to the Mach number 
Ms and the adiabatic index of a non-relativistic ideal gas 7 = 5/3. The columns '(?'tot)Mc' 
and '(T2/To)p^q' show the result of the simulation, i.e., the average of the last 10 of the 
30 iterations, along with the la standard deviation. For Ms = 0.5, the Hugoniot adiabat 
doesn't apply, because the flow is subsonic, and no shock should form, corresponding to 
^'tot = 1 and T2 = Ti (values marked with the asterisk '*'). 

It is clear that, within statistical errors, the simulation n Models 1-9 reproduced 
the Hugoniot adiabat results for nonrelativistic hydrodynamic shocks. Note that the physics 
that the simulation is based on at this point involves only the isotropic particle scattering, 
the Lorentz transformations between the reference frame of the flow and that of the shock, 
and the requirement that the calculated momentum, energy and mass fluxes be conserved. 
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Figure 3.11: Iterative estimation of compression ratio. 



Figure 13.111 illustrates how the prediction of the self-consistent compression ratio 
for Model 9 starts off with a trivial initial guess of rtot = 1.001, rises up to the value 
predicted by the Hugoniot adiabat in 15 iterations, and stays there. 
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Figure 3.12: Momentum and energy conservation illustration. 
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Figure [3.121 confirms the conservation of momentum and energy across tlie sliock 
in Model 9. The top panel shows the flow speed with the compression ratio rtot = 4.0, the 
middle panel shows the momentum flux measured in the units of the upstream flux, and the 
bottom panel shows the energy flux. These fluxes remain constant throughout the shock 
for Model 9 (the dashed thick line). 

Now let's analyze the results for Model 10, which is the same as Model 9, but with 
the shock transparent to the particles, so that the process of diffusive shock acceleration can 
occur. The Hugoniot adiabate is not applicable in this case, and the corresponding values in 
Table [3T] are marked with a double asterisk '**'. Indeed, the values of the compression ratio 
and the downstream temperature determined by the simulation and shown in Table 13.11 are 
^tot ^ 3.09, and T2/T0 ~ 230. These values are is lower than the prediction of the Hugoniot 
adiabat for a shock with the sonic Mach number Mg = 30. In Figure [3.111 the thick solid 
line demonstrates that the iterative procedure given by Equation (j3.58p has converged, but 
the fluxes shown with the solid lines in Figure 13.121 are not constant, meaning that the 
derived solution is not physical. This is the effect of particle acceleration studied by Ellison 
and co-workers using the Monte Carlo model, and the commonly accepted solution is that 
a shock precursor must form, i.e., the upstream flow speed u{x < 0) must decrease towards 
the shock. The procedure that models it is demonstrated in Section 13.1.81 

3.1.8 Nonlinear structure of the shock precursor 
Theory 

In the test-particle limit, the aftermath of particle acceleration by a shock with 
the compression ratio r is the power law spectrum of accelerated particles every point 
in space, f{x,p) oc p~^. The power law index of the spectrum for a strong shock with 
r = 4.0 is s = 4.0, which has the unphysical property that, if it were to stretch from p = 
to p = 00, the pressure (and the internal energy) of the particles with such a spectrum, 
P oc J pf{p)p^ dp a J p^^ dp, would logarithmically diverge at the high momentum end. In 
reality, of course, the spectrum is limited by the maximum momentum determined by the 
size of the shock or the acceleration time. However, the logarithmic divergence of pressure at 
high energies means that there may be a large amount of energy in the highest momentum 
particles. Of course, the actual amount of energy is determined by the rate of particle 
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injection. 

Initially, the shock in the simulation doesn't have a self-consistent structure be- 
cause it starts with an unmodified shock and ^p{x) is overestimated at all locations where 
accelerated particles are present (see, e.g., the plots for Model 10 in Figure [H?T^ . Therefore, 
it must choose u{x) to reduce the mismatch between the local momentum flux and the far 
upstream value of it <I>po for x < 0, as described by Equation (j3.42p . It can be done by 
calculating 

n'{x) = u{x)+C.^^^^^l^^, (3.59) 
Pouo 

where u'{x) is the predicted flow speed for the next iteration, and C is a small positive number 
(typically around 0.1), characterizing the pace of the iterative procedure. The value of the 
parameter ( is randomly chosen between and a finite value for reasons similar to those 
described in Section 13.1.71 If magnetic field amplification is invoked, then at this point the 
simulation also refines its estimate for the particle diffusion coefficient (see Section [3. 3p . 

The predicted u{x) and D{x,p) are then used in a new iteration where particles 
are injected and propagated. The calculated CR pressure, momentum flux, etc. are then 
used to refine the guesses for u{x) and D(x,p) for the next iteration, along with the guess 
for the compression ration r^Q^. This procedure is continued until all quantities converge. 

Test of implementation 

In order to test the implementation and the effects of the precursor smoothing, I 
complemented Model 10 with the procedure (j3.59p . I had to reduce the maximum value 
of C from 0.3 to 0.1 and make many more iterations in order to restrict the self-consistent 
compression ratio rtot- As Figure [3.131 shows, the iterative procedure converged to a value 
^'tot ~ 11, much higher than the Hugoniot adiabat predicteco- 

In Figure [3.141 1 show the spatial structure of the smoothed shock. In the top panel, 
a reduction of the flow speed from uq in the upstream region, x < 0, is apparent. Let us 
refer to the smoothed region as a shock precursor. A subshock at x = has a compression 
ratio, Tsub = u{x = —0)/u{x = +0) ~ 2.5. Momentum flux shown in the second panel 
is conserved within a few percent, and its deviation from conservation is statistical (the 



This and other nonlinear effects axe discussed in 1431. 
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Iteration 



Figure 3.13: Search for the self-consistent compression ratio. 



shown plots are an average of 14 iterations at the end of the 500 iterations leading to a 
consistent solution). The energy flux (third panel from top) drops at the upstream free 
escape boundary, x = —80 rgo, by 60%, which is explained by particle escape. Downstream 
of the free escape boundary, the energy flux is almost constant (within statistical deviations). 
The thin dashed line in this panel shows the self-consistent value of energy flux accounting 
for particle escape, as determined by equation ()3.53p : as one can see, the actual value of 
is in excellent agreement with this quantity. The bottom panel shows the constituents 
of the momentum flux <I>p, the dynamic pressure pu^, the thermal particle pressure Pth 
and the cosmic ray pressure Per- One can easily see that the shock is dominated by the 
accelerated particles in this case. 
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Figure 3.14: Precursor smoothing for momentum and energy conservation. 
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3.1.9 Summary 

I would like to conclude this lengthy section with a summary. The simulation of 
nonlinear shock acceleration presented here features a Monte Carlo code of particle transport 
and an iterative procedure for deriving a self-consistent shock structure. 

The tests I presented here confirm that: 

• Particle propagation simulated by the Monte Carlo method is diffusive, and that it 
reproduces qualitatively and quantitatively the expected behavior of such particle 
transport; 

• The simulation reproduces the well-known results for hydrodynamic shocks (i.e., the 
Hugoniot adiabat), if particle acceleration is artificially blocked; 

• Accelerated particle spectra predicted by the simulation in the test-particle regime 
agree with the solid predictions of analytic models of test-particle acceleration; 

• When the feedback of accelerated particles in the nonlinear regime is accounted for, 
the simulation obtains a stable solution that conforms with mass, momentum and 
energy conservation laws. 

The methods of the Monte Carlo simulation of particle transport and the idea 
of iterative derivation of nonlinear shock structure presented here are not original to this 
dissertation, they were developed by Ellison and co-workers. I will not elaborate any more 
on the subject of nonlinear shock acceleration and refer the reader to the literature for more 
information (e.g., [72j and [92j). 

The original part of the model - magnetic field amplification and self-consistent 
particle transport - will be presented in the rest of this Chapter. The reason I presented 
and tested the Monte Carlo part of the model in such detail is that the actual computer 
code used in the research was written by the author of this dissertation and had to be tested 
to confirm that it reproduces the well known, previously established results. Besides that, 
some details of the currently employed methods did not appear in our publications, and I 
would like to make a record of these details here. 
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3.2 Magnetic Field Amplification 

In this section I will describe two theoretical models of magnetic turbulence am- 
plification by streaming particles available in the modern scientific literature. Then I will 
proceed with a generalization of these models for the problem of nonlinear DSA, and present 
the analytic description of magnetic field amplification adopted for the model, and its nu- 
merical solution. I will conclude by describing the feedback of the amplified magnetic 
turbulence on the plasma flow. 

The physical conditions in which the instabilities take place are representative of 
the conditions in a collisionless shock precursor. A fully ionized plasma is moving at a speed 
u{x) from the far upstream, unshocked region (x — > — oo) towards the subshock located at 
X = 0, where it gets non-adiabatically compressed. We assume that a pre-existing uniform 
magnetic field Bq parallel to the plasma flow fills the space. The instabilities are induced by 
the accelerated particles produced by diffusive shock acceleration, and described by /(x, p). 
These particles are subject to diffusion in the plasma and advection along with it. Therefore, 
in the reference frame locally co-moving with the plasma, the CRs appear to move against 
the bulk flow, away from the subshock. The CR density and pressure increase from at 
x — > — oo to a finite value at x = 0, and so does the diffusive current of CR measured in the 
plasma reference frame ^. 

We describe the fluctuations of magnetic field in the plasma by the energy spectrum 

of turbulence, W{x, k). The latter is a quantity such that W{x, k)Ak is the volume density 

of turbulent energy (i.e., the energy of the waves, including the magnetic field energy and 

that of the associated stochastic plasma motions) in the waveband Ak. This means that, 

instead of the the spatial structure of the turbulent magnetic fields, the simulation only 

follows the evolution of its Fourier transform, locally calculated over a large enough volume, 

and averaged over a large enough time interval. This approach relieves our model from the 

computational expenditure of PIC plasma simulations, because the latter have to have a 

spatial resolution finer than the size of the smallest turbulent vortex, while our simulation 

gets away with a resolution that is more coarse than the largest turbulent harmonic by 

describing the processes on smaller scales statistically. Obviously, this advantage is gained at 

^Speaking of CR current and CRs streaming, I will always mean the apparent drift of CRs in the plasma 
reference frame. 
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the cost of having to rely on theoretical models of processes in plasmas instead of performing 
a numerical experiment based on more fundamental physical principles. 

Magnetic turbulence amplification is modeled under the assumption of a steady 
state situation, i.e., the time derivatives d/dt in the respective equations are set to 0. It 
is also postulated that in the interstellar medium (i.e., far upstream), there exists seed 
turbulence, which is expressed by the boundary condition 



W{-oo,k) = < 



47r In (/Cmax//Cmin) (3.60) 

0, otherwise. 



Expression (j3.60p describes a far upstream seed turbulence with a power law spectrum, 
W oc k~^, normalized so that Be{f[x — > — oo) = ABseed, where Ai?seed is a parameter of 
the model representing the assumed effective magnetic field of the seed turbulence. The 
values kmin and fcmax limiting the wavenumber range of the seed turbulence spectrum are 
also parameters of the model. Normally, the code will choose them so that for particles of 
all energies found in the simulation, the resonant wavenumber (defined later) is between 

^min and A^jnax- 

The effective magnetic field Bcsix) at any point is a quantity that I define as 



Stt 2 



W{x,k)dk. (3.61) 



The fraction 1/2 before the integral in the right-hand side of Equation (|3.6ip expresses 
the assumption that one half of the turbulence energy is contained in the magnetic field 
fiuctutaions, and the other half is carried by the stochastic fluctuations of the plasma velocity 
associated with the waves. The factor 1/2 is exactly correct when the turbulence is purely 
Alfvenic (see Section I3.2.9P . The expression (I3.6ip is therefore a good approximation if the 
turbulence is generated by the resonant streaming instability (see Section 13.2. ip , but does 
not apply, for example, to the waves generated by the nonresonant Bell's instability (see 
Section [322]) • For the latter case, the relationship between the plasma velocity fluctuations, 
6u, and the magnetic field fluctuations, 6B, can be inferred, for example, from Equation (17) 
in [To]. It can be shown that for wavelengths at the peak of the ampliflcation rate, k = 
kc/2 (deflned in Section I3.2.2p magnetic field contains 3/4, and velocity fluctuations - 
1/4 of the total energy density W{x,k). However, I use the '50/50' distribution of the 
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turbulent energy between magnetic and kinetic fluctuations, expressed by the factor 1/2 in 
Equation (j3.6ip . for turbulence produced by any source. I do so for the following reasons. 
First, considering the much larger uncertainties in some other factors of the model (e.g., 
the diffusion coefficient), an adjustment of the expression (j3.61|) to account for the nature 
of the turbulence would be a minor correction. Second, nonlinear turbulent processes like 
dissipation and cascading may change the energy distribution between magnetic and velocity 
fluctuations, and there exists no analytic description of this process adequate for the strong 
turbulence considered in this model. 



3.2.1 Resonant cosmic ray streaming instability 

Charged particles streaming along a uniform magnetic field are able to resonantly 
amplify Alfven waves traveling along the same field (|8i nil8| 11261 19|. etc.). The mechanism 
of Alfven wave amplification is similar to that used to amplify electromagnetic waves in the 
electronic device known as the traveling-wave tube [63j. Alfven waves exchange energy with 
fast particles of resonant momenta. If the particle distribution is anisotropic, Alfven waves 
traveling in one direction get amplified, and waves traveling in the opposite direction get 
dampened by the interactions with the particles close to the resonance. I will not describe 
this instability in detail and refer the reader to the sources cited above. 

Assuming a steady state and that the particle distribution is controlled by advec- 
tion of the fiow and resonant turbulent diffusion (see Section I3.3.2P , and that the generated 
waves are a weak perturbation of the uniform magnetic field, AB <^ Bq, the evolution of 
W{x, k) may be described (see |82j and references therein) as 

^^^^ dW{x^,k) ^ ^^^^^^^ ^^^^^^ 

if all other effects accompanying wave generation are ignored. Here the growth rate 



■n / 7\ dPcr{x,Pres) 

r^esix, k) = VA 

OX 

and 



dk 



^ (3.63) 



W{x,k)' 



^k = 1. (3.64) 

eBo 

defines the resonant momentum, pj-^s- Here va is the Alfven wave speed, and Pcr{x,p) is 
the spectrum of particle pressure, normalized so that PcrAp is the pressure of accelerated 
particles in the momentum range Ap. 
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It was usually assumed that the fluctuations grow until IS.B ~ Bq, after which 
the instability saturates (e.g., |93|)- Recently, Bell and Lucek ^11 j suggested that, if the 
instability can grow beyond this point, to AS ^ Sq, the observations of large magnetic 
fields in some SNRs can be explained by generation of magnetic fields in the process of DSA. 
Simulations done by the same authors [86] support such a possibility. Here, I adopt Bell 
and Lucek's idea and assume that the streaming of CRs is able to produce strong magnetic 
field fluctuations. 

As the perturbations grow and reach /S.B > Bq, however, it is likely that waves with 
wave vectors k not aligned with Bq will be generated, due to local CR pressure gradients 
along the total B = Bo + AB. With AB > Bq, it becomes impossible to predict the average 
value of the transverse pressure gradients and the resulting magnetic field structure without 
knowing the relative phases of different wave harmonics. The problem is further complicated 
by the fact that this longitudinal, compressible turbulence may produce a strong 2nd order 
Fermi particle acceleration effect which, in turn, can damp the longitudinal fluctuations 
(see, for example, |105| ). 

These complications place a precise description of plasma turbulence beyond cur- 
rent analytic capabilities. However, valuable conclusions about MFA in efficient DSA can 
be made by considering the two limiting cases of the resonant instability development in 
the nonlinear regime. The first assumes there is no longitudinal turbulence, in which case 
the wave growth rate is determined by the Alfven speed in the non-amplified field Bq. This 
gives a lower bound to the growth rate. The upper limit assumes that the turbulence is 
isotropic, in which case the growth rate is determined by the Alfven speed in the much larger 
amplified field i?efr (defined in Equation (13.6ip ). The real situation should lie between these 
two cases, and while I consider these limits, I do not explicitly include second-order Fermi 
acceleration in the calculations. Section 14.11 describes the parameterization of the growth 
rate of the resonant instability that encompasses the minimum and the maximum growth 
rate in the regime of strong fluctuations. 
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3.2.2 Bell's nonresonant instability 

A nonresonant instability, theoretically described by Bell in 2004 ([10]). occurs 
when a strong external electric current of CRs is put through the plasma. The instability 
develops because the thermal plasma must provide a current in response to the external 
current of streaming CRs, in order to maintain quasi-neutrality. This current makes certain 
MHD modes unstable; these modes can be described as driven circularly polarized Alfven 
waves. Again, I will not discuss the details of the instability; the reader may find more 
information in [10^ 1130^ Hj . 

In the linear regime (AB <^ Bq) the growth of the unstable modes can be described 
by the following equation: 

u{x)^^^^r^ = rnr(x, k)W{x, k). (3.65) 
ox 

The dispersion relation of waves subject to Bell's instability is 

^2_^2^2^MM=0, (3.66) 
cp 

where lo and k are the frequency and the wavenumber of the generated waves, and jd is 
the diffusive current of CRs directed along the magnetic field Bq. The frequency lo has an 
imaginary part when 

k<kc = (3.67) 
cpv\ 

and the reasoning leading to (I3.66P is applicable when the wavelengths of the generated 
waves are shorter than the smallest energetic particle gyroradius in the system, rgi: 

— <k. (3.68) 

rgi 

Along with its applicability conditions (j3.67p and (|3.68p . the growth rate of the energy of 
the waves Tarik) = 2Imu;(A;) is 



r 

-■- nr 



0, otherwise. 

Here va = Bq/ ^JAiip is the Alfven speed, and the critical wavenumber kc is 

k, = (3.70) 

cpv\ 
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This instabihty has been studied theoretically and using MHD and PIC simulations 
(p ^ 11301 l5l HI [TOi] ). which show that it is capable of generating large magnetic fields, and 
may even dominate the resonant CR streaming instability in young shocks of SNRs [99]. 

It is informative for future reasoning to point out the dependence of Fnr on k: 
the rate Fnr becomes non-zero at large wavelengths at A; = l/^gi, and then grows as 
towards the smaller wavelengths, until it peaks at k = kc/2, and then rapidly falls off down 
to zero at k = k^- 

3.2.3 Nonresonant long- wavelength instability 

Other instabilities possibly leading to magnetic field amplification may exist in 
an interstellar plasma in vicinity of a particle accelerating shock. Bykov and Toptygin 
|30j suggested a model in which streaming cosmic rays may amplify waves in plasma with 
wavelengths much larger than the gyroradii of the particles. 

The model presented in [30] requires the presence of a neutral component in the 
plasma (i.e., unionized hydrogen atoms) that suppresses the transverse conductivity. The 
Balmer series lines of neutral hydrogen have been observed in the emission spectra of for- 
ward shocks of Type la SNRs (e.g., in SN 1006 and Tycho's SNR [77]. See also the references 
in [60]). It is also possible that short-scale turbulence may act similarly to neutral plasma 
component at suppressing the transverse conductivity, which will make the model of nonres- 
onant long-wavelength instability applicable to fully ionized plasmas as well (A. M. Bykov, 
in private communication, and [27J, in press). 

I have not incorporated this effect into the simulation, because the details of this 
model are being developed. However, if this instability operates in the precursor of a 
collisionless shock, it may have a very strong impact on the maximum momentum of the 
accelerated particles due to its long-wavelength nature. 
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3.2.4 Evolution of turbulence in a nonlinear shock — other effects 
Flow compression 



const, and do not take into account the possibility of a changing flow speed in the precur- 
sor. In the geometric optics approximation, the propagation of waves in a medium with 
changing properties may lead to modulation of the wavelengths and of the amplitudes of 
the waves. Diff'erent instabilities generate different types of waves, which may evolve dif- 
ferently. Namely, the resonant CR streaming instability in the quasi-linear approximation 
AB <^ Bq generates Alfven waves propagating against the flow, and Bell's nonresonant in- 
stability generates almost purely growing harmonics that can be described as driven Alfven 
waves [lOJ . Considering that compression ratios in strong SNR may reach rtot ~ 5 — 15 (see, 
e.g, observational arguments [1251 [33] . and theoretical predictions [72l [92]), these effects 
may change the wavelengths and amplitudes of the generated waves by a large factor, and 
should be considered. 

1 propose to include the effects of flow compression in the model by adding the 
corresponding terms to the equation of turbulence evolution. Consider the equation 



Given a boundary condition W{xq, k) = W(){k), one can readily solve this equation for x > 
xq using the method of characteristics: see Equation (j3.104p . The ratio u{xq) / u{x) in this 
equation is the factor by which the plasma is compressed, because due to Equation (j3.4ip . 
p{x)/p{xo) = u{xq)/u{x). Equation (13.104p shows that the parameter a determines how 
the amplitudes of the waves react to compression, namely, W (x p". For Alfven waves, 
the correct value of this parameter is a = 1.5 (see the wave generation equations and 
the corresponding explanation in [122J). The parameter /? describes what happens to the 
wavenumber of the waves as they propagate in the compressing medium, that is, k (x p^. 

Dissipation 

The amplifled turbulence may be dissipated through collisional and/or collisionless 
mechanisms and these include: (i) linear and nonlinear Landau damping (e.g., pi 1801 [T20| 
I128j ). (ii) particle trapping (e.g., [94j), and (iii) ion- neutral wave damping (e.g, [391 130]). 



Equations of the form (j3.62p and (j3.65p generally apply to a uniform flow u{x) 




(3.71) 
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Existing analytic descriptions of MHD wave damping rely on the quasi-linear approxima- 
tion AB ^ Bq, which is inapplicable for strong turbulence, and numerical models with 
varying ranges of applicability have been proposed which offer a compromise between com- 
pleteness and speed (e.g., [TU l \TT2\ ll3Uj ). Because no consistent analytic description of 
magnetic turbulence generation with AB > Bq exists, and because an numerical (MHD or 
PIC) description of this process in the framework of non-linear DSA is very computationally 
expensive, we propose a parameterization of the turbulence damping rate. In doing this, 
we are pursuing two goals. First, we make some predictions connecting cosmic ray spec- 
tra, turbulent magnetic fields and plasma temperatures, which, in principle, can be tested 
against high resolution X-ray observations in order to estimate the heating of the thermal 
gas by turbulence dissipation. And second, once heating is included in our simulation in a 
parameterized fashion, we will be ready to implement more realistic models of turbulence 
generation and dissipation as they are developed. 

The heating of the precursor plasma by dissipation modifies the subshock Mach 
number (e.g., |44 1I123| ) and this in turn modifies injection. The overall acceleration efficiency 
and, of particular importance for X-ray observations, the temperature of the shocked plasma 
(e.g., [Ml EHl [5ll ) will depend on wave dissipation. 

The generalized way of including the dissipation of turbulence into the simulation 
is introducing a corresponding term into the equation of turbulence evolution: 

dW 

u— = G-L, (3.72) 

where G stands for the rate of growth of the instability driving the turbulence (i.e., G = 
TW), and L = L{x, k) is the rate of turbulence dissipation (measured in ergs-cm~^s~^). In 
the simulation, I allow for one of two prescriptions for the dissipation rate to be realized. 

In the absence of a better model, one may assume that the dissipation rate L is 
proportional to the amplification rate G, i.e. 

Lp = ohG, (3.73) 

where an is a number between and 1. Equation (j3.73p is a mere parameterization of 
the dissipation rate, in which an is the fraction of the instability generation rate that is 
assumed to go directly into particle heating rather than magnetic fields {in Lp, the subscript 
F stands for 'fraction'). In particular, an = corresponds to no dissipation, and an = 1 
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describes a situation where all the turbulence generated by a particle streaming instability 
immediately gets dissipated and transformed into heat. 

Another prescription for the rate of dissipation is 

Lv = ^k^W. (3.74) 

The wavenumber at which the dissipation begins to dominate, kd, is identified with the 
inverse of a thermal proton gyroradius: 

= , (3.75) 

c^/rripkBT 

where rup is the proton mass, ks is the Boltzmann constant and T = T[x) is the local 
gas temperature. This prescription is a relaxation-time approximation, defined in such a 
way that at k = kd, the dissipation time equals the period of an Alfven wave with the 
wavenumber k, and the /c-dependence is L oc k'^. The k'^ dependence of the dissipation rate 
is based on the assumption that viscosity (i.e., magnetic viscosity in this case) drives the 
dissipation (see, e.g., [120]). 

Influence of turbulence dissipation on thermal plasma heating 

Dissipation of turbulence acts as an energy sink, in which the magnetic and kinetic 
energy of turbulent fluctuations are transformed into the internal energy of the thermal 
particle gas. This means that, in order to conserve energy, the appearance of the term L in 
the equation of turbulence evolution must be accompanied by the corresponding correction 
to the equations of motion of the thermal plasma. The way to incorporate the thermal 
plasma heating due to turbulence dissipation was shown in [93], who derived the equation 
of thermal pressure evolution in the shock precursor: 

:i^^{P..p-)=L. (3.76, 

Here the ratio of specific heats of an ideal nonrelativistic gas is 7 = 5/3. For L = 0, equation 
(j3.76p reduces to the adiabatic heating law, Pth ~ p"^ and, for a non-zero L, it describes 
the heating of the thermal plasma in the shock precursor due to the dissipation of magnetic 
turbulence. The fluid description of heating given by equation (j3.76p . while it doesn't 
include details of individual particle scattering, can be used in the Monte Carlo simulation 
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to replace particle scattering and determine heating in the shock precursor. This merging 
of analytic and Monte Carlo techniques, or Analytic Precursor Approximation (APA), is 
described in detail in Appendix B of |123j . and briefly summarized below. 

When the heating rate, L, becomes available from the solution of the turbulence 
growth equation (j3.86p . the code solves (|3.76p and substitutes the solution, Ptii{x), for the 
thermal pressure calculated from particle trajectories. It is done in the upstream region 
up to the point xapa, at which thermal particles are subsequently introduced into the 
simulation for the next iteration. In order to include the effects of heating in the model, we 
must introduce thermal particles at xapa as if they were heated in the precursor, i.e., their 
temperature r(3;APA) must be determined by (j3.76p and the ideal gas law: 



kBno{uo/u{x))' 

The simulation therefore chooses the magnitude of every introduced particle's momentum 
p in the local plasma frame distributed according to Maxwell-Boltzmann distribution with 
temperature T determined by (I3.77P at x = xapa- As long as the local sonic Mach number 
at this location is large (i.e., Mgi > 3), it can be done using the results of the section [3.1.51 
for the angular distribution of the introduced particles. If the AIsi < 3 (which does not 
usually happen in a self consistent solution), then the results described in Appendix B may 
be applicable. 

Spectral energy transfer 

Observations of turbulence, including the MHD turbulence in the interplanetary 
plasma, often report spectra that look like power law functions of k over many decades. 
This phenomenon can successfully be explained by spectral energy transfer (cascading). 
After the energy has been generated by an instability on some dominant spatial scale, 
nonlinear motions in the turbulent fluid cause splitting and merging of the turbulent vortices 
(i.e., a cascade), leading to a re-distribution of energy between different scales. This way, 
turbulence initiated by large-scale vortex formation due to an external power source can 
cascade into smaller vortices, producing a power-law distribution of energy in the so-called 
inertial range (i.e., the interval in A;-space where the turbulence spectrum is populated by 
cascading rather than directly by the instability). This cascade continues until the size of 
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the vortices is small enough so that dissipation [e.g., (j3.74p ] terminates it by converting the 
energy of motion into heat in the so-called dissipative region of fc-space (e.g., |17]). 

There are various ways to describe spectral energy transfer (see, e.g., |96)). One 
of the simplest methods, listed in |96j as the Kovazhny hypothesis, involves a dimensional 
analysis argument. If one writes the equation of turbulence evolution in the inertial range 
as 

then by the physical meaning, 11 is the flux of energy through fc-space towards larger k. 
Assuming that 11 is a product of powers of the minimum set of relevant quantities, one can 
find the simplest form of the corresponding cascading rate. That is, if 11 = W°'k^p'^, then 
the only combination of a, b and c that gives 11 the correct units is 

As we will see later, this cascading rate gives a stationary solution W oc k~^^^, which is 
known as the Kolmogorov spectrum, and the corresponding cascade will be referred to as 
Kolmogorov-type cascade (here denoted by the subscript 'K'). 

When MHD turbulence is considered, the simple dimensional argument shown 
above does not work because the magnetic field is another relevant quantity. There are two 
approaches to describing nonlinear effects (spectral energy transfer) in MHD turbulence. 
One, proposed by Iroshnikov and, independently, by Kraichnan [70l[78], treats the MHD 
turbulence as weakly interacting plasma waves that can undergo mergers and splitting. 
The bottom line of this approach is that a stationary spectrum W oc k~^^'^ is predicted. 
Because 5/3 and 3/2 are so close, it is difficult to distinguish between the two indices in the 
analysis of observations. Goldreich and Sridhar |64j point out that the MHD turbulence 
is inherently anisotropic (even if there is no mean magnetic field, the effective field of the 
large scale harmonics can play its role for the processes in the inertial interval), and the 
weak-turbulence approach is not applicable. These authors proposed another theoretical 
approach: they suggested a certain anisotropic damping rate and postulated a critical 
balanced state, which allowed them to derive an anisotropic turbulence spectrum. Their 
results predict that harmonics with wavevectors transverse to the uniform magnetic field 
experience a Kolmogorov-like cascade, while the cascade in wavevectors parallel to the field 
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is suppressed. The waves generated with streaming instabihties are transverse; therefore the 
diffusion coefficient for particle transport parallel to the flow depends on the wavenumbers 
parallel to the magnetic field. Biskamp [T7| shows that the Goldreich-Sridhar spectrum for 
parallel wavenumbers may be expressed as oc A;|| ^^"^ . 

We can find the corresponding cascading rate, such that 11 = W"" which 
would lead to a steady state spectrum with W oc k~^/'^ . From the dimensional argument, 

Hgs = T^2/3^5/3^1/3^5/3_ ^gg^^ 

One may do a simple estimate and compare the Kolmogorov and the Goldreich-Sridhar 
cascading rates: 

IIgs _ P72/3fc5/3pl/3^5/3 ^ / ^2 X 5/6 
~ p^3/2^5/2^-l/2 - \A-KkW J ■ 

Transition to turbulence 



(3.81) 



One may pose a relevant question: at what point do the linear plasma waves 
acquire the nonlinear behavior that leads to their cascade and dissipation at short wave- 
lengths? We assume that it happens when some of the waves reach strong amplitudes, i.e., 
AB(fc) ^ Bq. In terms of the quantities that we use to describe the turbulence spectrum, 
I postulate that if, at the coordinate x, there is a wavenumber k such that 

hw{x,k)>^, (3.82) 

then downstream of this coordinate, the turbulent cascade and dissipation start (i.e., a 
transition to turbulence occurs). In accordance with that, upstream of this coordinate, the 
energy transport 11 and the dissipation rate L are both set to zero. 

Anisotropy relaxation 

The resonant streaming instability of Alfven waves amplifies the waves traveling in 
the direction of the diffusive particle stream (i.e., in the upstream direction) and damps the 
waves traveling in the opposite direction. The distribution of energy between the upstream 
and downstream traveling waves may be important for some applications: for example, the 
mean speed of scattering centers, if it is not negligible compared to the flow speed, calls 
for the appropriate reference frame transformations for particle scattering. This may be 
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important for low Alfven Mach number shocks. In the strong, fast shocks of young SNRs, the 
generated waves predominantly travel upstream, but for older shocks, nonlinear interactions 
between upstream and downstream traveling structures may lead to the appearance of 
downstream traveling waves (e.g., [11]). 

Following [11], one can notionally separate the turbulence spectrum as in the 
plasma frame 

W{x,k) = U-{x,k) + U+{x,k), (3.83) 

where U- is the spectral energy density of waves traveling upstream, and [/+ - that of 
the downstream-directed waves. The equation of turbulence growth due to a streaming 
instability, accounting only for the wave advection, growth and the nonlinear interactions 
between the waves traveling in different directions can be written (see also |122) ) as 

[u{x) - va]-^U. = r,esU- - VAk iU- - U+) ; (3.84) 
[u{x) + va]-^U+ = -r,esU+ + VAk ([/- - U+) , (3.85) 

where = cp/{eBQ) and va is the Alfven speed. The factor uzizVA represents the fact 
that the considered waves travel at a velocity va with respect to the plasma along the 
magnetic field. The term proportional to C/_ — C/+ in the right-hand side describes nonlinear 
interactions between the oppositely directed waves that lead to isotropization of the wave 
spectrum (i.e., to C/_ = [/+) with a relaxation time of about the Alfven wave period. This 
effect may be important for weaker shocks. 
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3.2.5 Generalized model of magnetic turbulence amplification 

Considering the effects described above (except for tfie interactions with the waves 
traveUng downstream), let us write the equation of turbulence spectrum evolution in the 
following parameterized form: 

+ - ^'i h'i) - -"^ = »• 

and assume that a boundary condition is given at the coordinate x = xq in the form 

W{xo,k) = Wo{k). (3.87) 

The coordinate xq is typically located far upstream of the shock, and the function VFo(fe) 
describes the seed turbulence spectrum that, we assume, exists in the unshocked interstellar 
medium. 

In Equation p.86p . the first term describes the advection of turbulence with the 
flow. In the Lagrangian view, one may think of the turbulence amplification process as the 
evolution of a matter element advected towards and across the subshock, compressed and 
penetrated by cosmic ray flux on the way, which leads to a buildup of stochastic magnetic 
fields in this element. In the Eulerian perspective, this term represents the full derivative 
of W with respect to time, d/dt = d/dt + ud/dx, with the local derivative d/dt set to zero 
to model the steady-state solution. 

The second term, proportional to Og, represents the effect of plasma compression 
on the amplitude of the waves, as described in the previous section. The parameter Og 
measures the degree of this effect. With all other terms set to zero, equation (j3.86p has 
the solution W oc i.e., the amplitude of the waves grows proportionally to the power 

ag/2 of the plasma density. For Alfven waves, a = 1.5. 

The third term, containing Pg, describes the effect of compression on the wavenum- 
ber of the waves. With all other effects inactive, (|3.86|) has the solution W{k) oc W{ku^)u^ , 
which means that the spectrum W{k) shifts in log k space, while preserving the normal- 
ization: J W{k)dk = const. Setting /3 to may be used to 'turn off' this effect in the 
model. 

The term that contains jg is the driving term of instability growth. The function 
G is G = TW, where the growth rate F can take on values: Fj-es from (|3.63p . or Fnr defined 
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by (|3.69|) . or the sum of the two, depending on which instabihty one wishes to consider 
in the model. The value 7^ = 1 can be used to 'turn on', and 7^ = - to 'turn off' the 
turbulence amplification for purposes of testing the code or making predictions relevant to 
the physics of shock acceleration. 

The parameter 5g in the fourth term of (j3.86p controls the rate of cascading. For 
the energy flux 11 = Hk given by (I3.79p . the quantity 6g is essentially the Kolmogorov 
constant, a factor that complements the dimensional analysis leading to the derivation of 
(j3.79p . and that should be taken from experiments or numerical simulations. There seems to 
be a universal value of the Kolmogorov constant (see [lllj for a review of experiments (note 
that this article has a different definition of the constant) and [65] for simulation results), 
6g = 1.6 — 1.7. As was mentioned earlier, MHD turbulence may have cascade properties 
different from those of hydrodynamic turbulence, and 11 may assume different forms, for 
example IIgs from (j3.80p . For lack of better knowledge, I will use the value 5^ = 1 to 
include turbulent cascading and 5g = to omit it from the model. 

Dissipation of turbulence is controlled by the last term in (j3.86p , and the parameter 
Eg can be set to 1 or to include or omit the dissipation. The function L can assume the 
parameterized form Lp from (j3.73p or in the form of viscous dissipation Ly defined in 

(|mp . 
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3.2.6 Analytic solutions for turbulence spectrum 

The expression (|3.86p is a nonlinear partial differential equation of first order for 
a function of two variables, W{x,k). Its particular solution is determined by the initial 
conditions (j3.60p and by the nonlinear dynamics of the system that couples W{x, k) to 
particle propagation, and the latter to the driving term G in (I3.86P and to the flow speed 
u{x) determined with the iterative procedures (j3.58p and (j3.59p . Therefore, the solver of 
p.86p has to be run after every Monte Carlo iteration to advance the solution towards 
self-consistency. 

A powerful tool for tackling first order nonlinear equations is the method of char- 
acteristics [85] . In fact, in some simple cases, i.e. when the terms in ()3.86p assume a simple 
form, analytic solution is possible. Although these simple cases are not directly applicable 
to the physical system we are studying, I would like to derive these solutions below, because 
not only do they reveal the influence of various effects on the solution, but they will also be 
used for testing of the numerical solver. 

Parametric form 

In order to apply the method of characteristics to Equation (j3.86p . let us re- write 
it, collecting the terms containing the partial derivatives of W, and assuming that 11 is 
given by ([3:79]) : 

dW ( n,du ?, W^l'^k^l'^\ dW 
+ [-^^''Tx + 2^^^V^ ) -dk = 

= iPa - - -6g ^^1^ + 7,G - e,L. (3.88) 

In the spirit of the method of characteristics, this equation can be written in a parametric 
form [85], describing x^ k and W as functions of a new parameter t: 

dx 

dk ^ ,du 3^ W^/^k^/^ 



dW _ 5 . 1^3/2/^3/2 

P 



IT - ~ "^^d^^ ~ 2^^ ni/2 +^9G-e,L (3.91) 
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The other two parametric equations form a system of nonhnear mutually dependent equa- 
tions, which cannot be solved in a closed form, but the solution can be expressed in a 
form suitable for analysis. Consider the substitution p[t) = W{x{t),k{t))u'^s-I^9(^x{t)), 
q{t) = k{t)u^3 (^x{t)) , and notice that {d/dx) = u~^{x{t)){d/dt). Then the above system of 
equations can be written as 



dx 
'dt 



: u, (3.92) 

dg _ 3 

dt ~ 2^' pV2 



(3-93) 



^ = _^5^^!^!^^-./2"/5. + (^^G-e,L)n"«-^«. (3.94) 
at I p^i'^ 

Equations ()3.92p , ()3.93p and (I3.94p are the desired parametric form of the generalized equa- 
tion of turbulence evolution (j3.86p . This form will be used in the numerical solver. The 
physical meaning of the parameter t is obvious from equation (13.921) : it is the time elapsed 
since a particular harmonic at A: = A;o started evolving at x = xq (corresponding to t = 0). 

I must point out that this system is strongly nonlinear, because the quantities G 
and L depend not only on the coordinate x and the wavenumber A;, but also on the values 
of "p and and not only locally, but also on the integrals of W oc p with respect to x and 
(7, via the transport and particle acceleration properties of the turbulence. 

However, assuming simplified expressions for G and L, we may obtain analytic 
solutions, as shown in the next two sections. 

Solution without cascades 

In the absence of cascading (bg = 0), equations (|3.92p . (j3.93p and (j3.94p are not 
coupled, and have the obvious solution 

/dx' 

xo 

q{t) = qo, (3.96) 

i 

Pit) = P0 + j ilgG-egL)u'^^'P^ dt'. (3.97) 
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or, in terms of k and W , 



dx' 
u{x') ' 



k{t)[u{x{t))f' = kMxo)f\ 
W{t)[u{x{t))r'~''' = W{xo,k^)[u{xo)r~^^ + 



(3.98) 

(3.99) 
(3.100) 



+ j {jgG{x',k')-egL{x',k')) [u{x')p~''' dt', (3.101) 

to 

where x' = x{t'), k' = k{t'). Finally, for the spectrum in terms of the original variables, 
W{x, k), we may write 



W{x,k) = Wo{k 



u{x) 

u{xq) 



u{x) 



1 o^g-Pg 



+ 



+ / {lgG[x\k 



u{x) 
u{x') 



u[x) 
u{x') 



ll^g 



u{x') 
u{x) 



ag-l3g 



dx' 



U[X' 



- (3.102) 



XQ 



Expression (|3.102p is the solution of equation (j3.86p with the boundary condition (|3.87p for 
(5g = (no cascading). 

In particular, setting = /?g = in ()3.102p . we get the solution describing the 
turbulence evolution with only amplification G and dissipation L accounted for: 



W{x, k) = Wo{k) + / {igG{x\ k) - egL{x\ k)] 



dx' 
u{x') 



Assuming the opposite, Ig = = 0, but 7^ and fig 7^ 0, one obtains 



(3.103) 



W{x, k) = WQ\k 



uix 



n(xo)) 



Pg' 



u{x) 



Oig-Pg 



(3.104) 



which describes the effect of compression on the plasma turbulence: the energy density 
increases in proportion to and the wavenumber grows as p^f (see also Section [3.2.4p . 



Integral form with cascades 

Now let us return to the parametric form of the turbulence evolution equation 
given by (|3.92p . (j3.93p and (|3.94p . This time let us not set 6g = 0, but try to derive a 
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solution that accounts for cascading. The last two equations are coupled via the cascading 
terms, so in order to get an integral form of the solution, let us express p as a function of q 
by dividing equation (|3.94p by equation (|3.93p . which is a correct operation because 5g ^ 0. 
This leads to: 

(3.105) 



dp__5p 2 1 _ 3„^/2 

dq~ 3 g 3^,1/2^5/2 5/^^^ ^^^^^ ' 



or 



d 
dq 



pq 



5/3\ 



3/2' 



which gives the dependence of p on g in the following form: 



-5/3 



3/2 1 



+ 



p'^/\-fgG{x' , k') - egL{x\ k'))u^"/\x')dq' 



-'9 J go 



2/3 



(3.106) 



(3.107) 



where x' = x{t'), k' = k{t'), and t' is the moment in time corresponding to q{t') = q'. 
Constants po and qo define the characteristic curve by its initial conditions as: po = 
W{xo,ko)u°'-f^{xo), and qo = kouf^ixo). Substituting (|3.1U7|) into ([3:93]) gives: 



d 

It L 



-2/3 



{vrnT)"'^ + j f p^(x')(7<,G(x', k') - SgLix', k')y/\x')dq' 



nl/3 



-'g J go 



(3.108) 



Expressions (|3.92p , ()3.107p and ()3.108p are almost a solution to equation (j3.86p . We collect 
them below: 



d_ 

dl 



dx 

-2/3 



.-a,/2-f^ -1/2 



5/3\3/2 



PoQo 



+ 



+ - 



p5(x')(7gG(x , k') - egL{x', k'))u^"/\x')dq' 



^9 J go 



1/3 



(3.109) 



(3.110) 



-5/3 



3/2 



+ 



1 r p'/\jgG{x',k')-egL{x',k'))u^^/\x')dq' 



■'a J go 



2/3 



(3.111) 



The system ()3.109p . (|3.110p and (|3.11ip . the integral (integro-differential) form of the sO' 
lution to equation (j3.86p . is rather unsightly for a physicist, so we will simplify it by set' 

ting Ug = Pg = Eg = 0, 'J g 



6g = 1, and assuming that u{x) = uq, p{x) = po, and 
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G{x,k) = Go(5d(A; — kc), where do is the Dirac delta function. This corresponds to the 
case where energy is supplied to the turbulence at a wavenumber kc throughout the spa- 
tial extent of the system. In addition, let us assume that there is no seed turbulence, i.e., 
Wo{k) = 0. Then (|mU]) and ([mi]) yield: 

-3/2 



k, 



k, 



W{x{t),k{t)) 



-2/3 



-2/3 



-2/3 



t < tc, 



-3/2 



(3.112) 



t > U 



.-2/3 



5/2 „2/3 

PoF^'\ t < t. 

^ X 1/3 15/2 

t 



po J 



Po{Fo + ^Y'. 
Po J 



(3.113) 



t>tr,. 



Here 



and 



WjxQ, fco )fco 
Po 



5/3 \ 3/2 



-2/3 



-2/3 



1/3 



(3.114) 



(3.115) 



Assuming Gq/pq ^ Fq (that is, the generation of the turbulence at the wavenumber kc 
overpowers the seed turbulence), the solution for t > tc (in other words, for k > kc) is: 

-3/2 

(3.116) 



k{t) 
Wix{t),k{t)) 



k, 



-2/3 



-2/3 



Po 

Gq 

po 



1/3 



5/2 



Po 



= k-'/Ht)Po (^X%-n7) 

Po J \Po J 



The expression for W{x,k) does not depend on k^, xq or W{xQ^kQ), and therefore it de- 
scribes explicitly the turbulence spectrum at A; > k^. Namely, shortward of k^ the effect of 
cascading leads to a formation of a power-law spectrum of turbulence W{k) oc A;^^/^, which 
is the Kolmogorov spectrum, as discussed in Section 13.2.41 This result can directly be used 
for testing of the numerical routine solving the equation (j3.86p . 
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3.2.7 Development of the numerical integrator 

In order to calculate the spectrum of MHD turbulence produced by the instabilities 
of the precursor plasma in the presence of the accelerated particle stream, the model solves 
equation p.86p . The driving term in this equation, G, is calculated using the information 
about particle streaming simulated in the Monte Carlo transport module. The numerical 
procedure that will be run in the simulation must solve Equation (j3.86p with arbitrary 
driving term G and with or without all the other terms in this equation, parameterized by 
ag, Pg, jg, 5g aud £g. I havc developed such an integrator, and the algorithm of integration 
is presented in this section. 

In brief, equation (j3.86p is solved by integrating the system of coupled first-order 
ordinary differential equations: (j3.92p . (j3.93p and (|3.94p . This system is derived using the 
method of characteristics, and its solutions for different values of ko are the characteristic 
curves. The numerical method used for integration is a finite differencing scheme (based 
on the implicit Gauss's method), with an adaptive step size in x-space and adaptive mesh 
refinement in /c-space. The implicit nature of Gauss's method is beneficial for the stability 
of the results, and is achieved with an iterative procedure. 

Here is the outline of the procedure. Integrating from x = — oo to x > 0, the 
scheme will make steps, being the number of grid planes. For every k-h'm, every 
spatial step from X(^i_i-^ to Xj-j) will consist of Ngui, substeps, enumerated by the index /, in 
which the code will propagate k and W from X(j_i) to the size of each substep will 
be adaptively chosen to ensure the stability of Gauss's method. After all A;-bins have been 
propagated from 2;(j_i) to the program will use the /c-grid modified by compression and 
cascading to project the amplified W onto the fixed /c-grid of the simulation at and 
then proceed with the step to the next grid plane. If the code finds that the evolved /c-grid 
has too large a spacing between some nodes, it will refine the problematic regions of the 
A;-grid at X(j_i) and repeat the integration of the equation. The scheme will keep refining 
the fc-grid at X(j_i) until the resulting /c-grid at X(j) is satisfactory (i.e., fine enough). 

Notation for this section 

An index in round parentheses, as in X(j), enumerates the x-grid plane, and one 
in square parentheses, as in ky^, indicates the number of /c-space bin. Subscripts without 
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parentheses (e.g., in pi) mean the number of the substep between and and 

superscripts in parentheses (e.g., q^"^^) are reserved for the number of the cycle in the 
iteration used to achieve the imphcitness of the method. 

For brevity, we re-write equations (|3.92|) . (|3.93|) and (j3.94p as 



dx 

It 
dq 

'dt 



u, (3.118) 



where 



Making a substep 



: 6gCq, (3.119) 
: -^6gCp+i-fgG-egL)u''^-^^, (3.120) 

3 pl/2g3/2 



The substeps wih be enumerated by the index /, so that qi and pi are the quantities 
q and p at the end of the l-th substep. The code starts making the / substeps by initiahzing 
the following quantities: 

xo = (3.122) 
to = 0, (3.123) 

qo = A:y](x(i_i))n^«(2;(i_i)), (3.124) 
po = VFy](x(,„i))u"«-^^(x(,„i)). (3.125) 

To make the l-th substep, let us first assign the following quantities: 

xi = xi-i + Axi, (3.126) 
ui = u{xi), (3.127) 
PI = p{xi). (3.128) 

The step width Axi will be initially (for I = 1) set as 

Axi = X(i) - (3.129) 

and if this attempted substep succeeds, there will be only one substep (/ = 1), after which 
the scheme will move on to the next grid plane i. If the scheme finds this substep too large. 
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it will choose a smaller substep. For the subsequent substeps we will set 



Axz = Xi ■ Ax/_i, 



(3.130) 



where Xi is a number, greater or smaller than 1, depending on whether the previous substep 
was estimated as too short or too long, as discussed later. The program can integrate (|3.118p 
to get: 

Axi 



tl 



ui 



(3.131) 
(3.132) 



To derive qi from qi^i and pi from pi-i, the code will need to use an iterative procedure 
in order to implement an implicit finite differencing scheme for solving ()3.119p and ()3.120p . 
The superscript (m) will denote the cycle of iteration, and it will run from as far as it 
takes for convergence, restarting as m = with each new /. The initial step in this iteration 
will be 



1i 



(0) 



Pi 



(0) 



qi-i, 
pi-i- 



(3.133) 
(3.134) 



and the subsequent iterations will be derived from 



1i 



(m) 



qi-1 exp 



Pi 



(m) 



iP 



7-1 




(m-1) 



Ati + 



(3.135) 
(3.136) 



The values of the above mentioned derivatives and of the quantity are discussed later. 
Before making the (m)-th iteration, the code must check whether the substep size Axi was 
small enough. It does so by comparing the arguments of the above mentioned exponentials 
to a pre-set number r/. The value r] = 0.01 seems to work well as the target step size. If at 
any step the arguments of the exponentials are greater than r], the (m) iteration terminates, 
the code chooses a proportionally lower Axi by setting Xi < 1, and tries making the l-th 
substep again. If the value of the arguments of the exponentials in (j3.135p and (j3.136p are 
by a factor of a few smaller than rj in all (m) iterations, then for the {I + l)-th substep the 
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code chooses Xi^i > 1 in order to speed up the integration. Choosing the spatial step size 
this way makes the scheme adaptive in x-space. 

The value p^, is used to tend to a nasty property of our equations: the cascading 
and dissipation terms eventually drive the solution to p{t oo) — > 0, which happens to 
be the boundary of the range of definition of some of the functions in the equations. For 
the analytic solution, it is not a problem, because at p = processes further decreasing p 
(cascading and dissipation) naturally cease. But in a numerical solution, there is a danger 
of marginally running into the p < region, if p is evolved with a finite differencing scheme, 
which will cause an error, because the factor p^/^ in some of the functions is not defined 
for a negative p. I eliminate the possibility of getting p < by evolving Inp instead of p 
with the finite differencing method. However, when p ^ 0, the program risks dividing by 
zero. To avoid zero values of p, I re-define the point at which the processes decreasing p 
stop: from p = to p = p-^. The solution in each bin subject to dissipation then converges 
to p = p^ instead of p = 0. The value of is chosen small enough so that it doesn't affect 
the physical solution, but large enough to be treated numerically without problems. 

One danger possible with an iteration on q^"^^ and like (j3.135p and (j3.136p 
is that the solution may find an attractor cycle around the equilibrium point instead of 
converging to it, in which case we may find ourselves stuck in an infinite cycle (the equi- 
librium point is the point at which q^^^ = ^'^ and pj™'' = p[™ ^^). Theory suggests 
that there is a finite domain of attraction around the attracting equilibria of this system, so 
all we have to do to ensure convergence in the end is perturb the solution occasionally. If 
the iteration is stuck in an attractor cycle, with the perturbation it usually jumps into the 
domain of attraction of the equilibrium point and converges (or finds another cycle, which it 
will be driven out of with a later perturbation). In practice, the code perturbs the solution 
whenever m equals a multiple of a large integer, for example, 1000. Then it adjusts q^"^^ 
and pf^^ only half way from q^"^ and pf^ to what (j3.135p and (|3.136p suggest (this 
going half way is the perturbation). Experience shows that this procedure successfully finds 
the equilibrium points of the above system of equations, thus yielding the implicit Gauss's 
integration scheme. 

The iteration deriving from and p^"^^ from p[™ will continue until 

it converges, that is, the relative difference between the values obtained at the previous and 



92 



the current step becomes small enough. Suppose it happens at step m = N^- Then the 
code will assign 



11 = Qi 



Pi = Pi 



iN,„) 



(3.137) 
(3.138) 



and increment I. As soon as the last substep is completed (x; = X(j)), the program names 
Nsub = I and assigns 

Qitfin) = QNsnt^ (3.139) 
Pit fin) = PN^^,- (3.140) 

Having q{tfin) and p{tfin) allows one to revert back to the physical quantities and assign 

ky]{Xi) = q{tfin)u{xi)-l^\ (3.141) 
W^^{x.i) = pitfinHxi)-"^^''^. (3.142) 



Calculating the derivatives 



In equations (I3.135P and (j3.136p . the derivatives of Inq and Inp are calculated, 
according to (|3.119p and ()3.120p . as: 



where 
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(m-l) 


dt 
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(m-l) 



Vg,i 



dPcr 

dx 



dp 



dk 



L 



(m-l) 



(m-l) \ -ag+fig 

Pi -pi,] Ui 



{xi,hk) 
H{xij,k) 



(3.143) 
(3.144) 



(3.145) 



(3.146) 



(3.147) 



' TD{Xi^uk) 

More details on evaluating quantities from ()3.146p and ()3.147p are given in the next sub- 
section. 
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Details of the the growth and damping rate calculations 

In expressions (|3.146p and (j3.147p the following notation is used: 



k 

Vga 



k 



(m-l) 



(m-l) -/3, 



Qi (3.148) 

The instability growth term, G, in the resonant case is determined by the gradient 
of CR pressure at the resonant momentum. The quantity Per is the pressure per unit interval 
of particle momentum, thus the factor \dp/dk\ in (I3.146p . To calculate the pressure gradient 
in such a way that the discontinuity of in p-space doesn't lead to a discontinuity of G 
in /c-space, 1 chose to average the pressure over a finite wavenumber interval A/c. The code 
sets A/c = 0.05A; and defines 

1 



k^''f\k) = k - ^Ak, (3.150) 
k"sf^\k) = k + ^Ak, (3.151) 

after which it can calculate the corresponding range of the particle momenta that interact 
with the current bin: 

eBo 



p"-'''\k) 
p^°'"{k) 



cU^ft{ky 

eBo 



(3.152) 
(3.153) 



ck"ght(^f^y 

Then the instantaneous gradient of the CR pressure that powers the instability (to the 
best of one's knowledge at the (m)-th iteration of the l-th substep from to Xj) can be 
estimated as 



dPcr 

dx 



dp 
dk 



1 



Ax 



{xi, k) 



P{{) (k) 



dp 



dk 



dp 



dk 



(3.154) 



where 



P(^-l) (k) 



P(:i) (k) 



dp 



dk 



dp 



Ak 



Pcv{x(i-l),p)dp, 



dk 



Ak 



Pcr{x{i),p)dp. 



(3.155) 



(3.156) 
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This allows us to calculate (|3.146p . Note that I used the grid nodes and as 

reference points for calculating the gradient. That is done because the CR pressure is 
evaluated in the Monte Carlo simulation directly at these locations. 
In the dissipation rate L, calculated in (|3.147|) . 

TD{xi,k) = -—, (3.157) 
k,{xi) = (3.159) 

Adaptive fc-grid 

I use the parametric form of the turbulence growth equation, in which the /c-grid 
evolves in time. It may happen that two fc-grid nodes that were adjacent far upstream 
will move apart significantly by the time the turbulence advects downstream. In practice, 
starting off at x = — oo with 80 A;-grid nodes equally spaced in log k space and spanning 
10 orders of magnitude of k, we are likely to get two adjacent nodes that move apart by 
several orders of magnitude (!) at x > 0. This makes it problematic to interpolate the wave 
spectrum W{x, k) between these two nodes, and, in fact, a lot of information about this 
fe-region is missing from the solution. An attempt to boost the /c-resolution by increasing 
the density of fc-grid nodes uniformly throughout fc-space leads to a significant increase of 
computation time and to the need to have tens of thousands of /c-nodes. 

To solve this problem, I use an iterative approach to the refinement of the k- 
grid. After the system is integrated to the code evaluates the A;-grid at this final 
point. If it finds two A;-nodes that are too far apart at x^j^ (by too far apart I usually 
mean AlnA; = lia(ky^ / ky > 0.5), it inserts a number of new nodes into the integration 
grid at and interpolates the seed turbulence spectrum into these nodes to repeat the 

calculation. Several (less than 10) iterations like that allow to get enough resolution in k- 
space throughout the system with minimal time (in practice, the whole computation takes 
a few seconds) and minimal memory (I usually have to have only a few hundred fc-bins). 
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Turning over in /c-space 

Equation (j3.90p . for (3g = 0, shows that cascading leads to the motion of a harmonic 
with wavenumber A; at a speed of Vk = l.bW^^'^k^^'^ p~^^'^ (for 6g = 1) in A;-space. This 
dispersion relation has an interesting feature: if the spectrum W{k) has a power-law shape, 
W oc k'^, then Vk is an increasing function of k for s > — 5, but a decreasing function of k 
for s < — 5. That is, for the parts of the spectrum in which it rapidly drops off with k (more 
quickly than k^^), the lower k harmonics increase their k faster than the greater k. In 
this situation the fast-moving low-fe harmonics may catch up and overrun the slow-moving 
high- A; harmonics. 

This situation is common for waves in gases and fluids, where the phase speed of 
waves increases with density or wave height. It leads to waves turning over in water, and 
to shocks in fluid dynamics, when viscosity is accounted for. Obviously, in this model the 
turning over of waves in A;-space doesn't have a physical meaning and simply reflects the 
limited applicability of the Kolmogorov cascade to steep wave spectra. However, straight- 
forward application of this model to non-linear particle accelerating shocks does lead to the 
turning over in A;-space in the numerical solution. 

The place where turning over is most likely to occur is the dissipative region of 
the spectrum. There the turbulence dissipation term, L, makes the wave spectrum drop 
off exponentially, creating the situation in which turning over is likely to happen. Another 
possibility is turning over in the inertial region, if the generation of waves occurs on top of 
a 'seed' spectrum, and the generated waves cascade faster than the seed waves. 

I ignore the wave turnover in /c-space in the dissipative region, assuming that it will 
not affect the energetics of the process very much. As for the inertial spectrum, the physical 
solution for a steady-state nonlinearly modified shock must not have wave turnover there, 
if the model is self consistent (otherwise we must conclude that the Kolmogorov cascade is 
not a good approximation for the plasma physics of self-generated turbulence). There is a 
natural property of the accelerated particle distribution in shock precursors that seems to 
help the situation, if resonant amplification of waves is assumed. Very far upstream, only 
the highest energy particles resonantly generate the smallest k waves. These waves start to 
cascade and would outrun the higher k seed waves, but as the plasma advances toward the 
subshock, it encounters lower energy particles, whose pressure builds up exponentially with 
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time. These lower energy particles should energize the higher k waves, facilitating their 
escape from the lower k waves pre- amplified farther upstream. This way wave turnover in 
/c-space may be avoided naturally due to the properties of particle accelerating shocks. 

3.2.8 Tests of the numerical integrator 

In this section I will present the tests of the integrator which compare the results 
of the numerical solution to the analytic solutions described above. All these test involve 
introducing a seed turbulence spectrum upstream, at x = xq < 0, and numerically integrat- 
ing Equation (j3.86p from x = xq to x = 0. In order to test and understand the effects of 
different processes parameterized by ag, (3g, 7^, 5g and £g, I executed several runs, in which 
some of these parameters were set to finite values, while the other were set to zero. 

First, I tested the effects of the compression of the flow: the increase in the ampli- 
tude and the wavenumber of the harmonics. At x = xq I introduced a Bohm seed spectrum 
with a Gaussian feature on top of it, located at A; = 10~^ r~Q (see the thin line in Fig- 
ure I3.15|) . I imposed a flow speed that drops by a factor of r = 10^ from x = xq to x = 0. 
Then the code solved Equation (j3.86p using Og = 1.0 and I3g = 2.0 (these values were used 
just for testing; physically justified values are discussed in Section r3.2.4p . The resulting 
spectrum at x = 0, shown with the thick line in Figure [3.151 agrees with one's expectation 
based on the analytic solution (j3.104p : the feature moved to the right, towards greater k 
by a factor of r^^ = 10^ and upward, to greater amplitudes, by a factor of r^f = 10^. Note 
that in the plots, the spectrum W{x, k) is multiplied by k, so a horizontal line represents 
the seed spectrum, W oc k'^ . 

The second test, illustrated in Figure [3.161 confirms that the amplification term 
proportional to 7^ in (13.860 is handled correctly by the numerical solver. The introduced 
seed spectrum (shown with the thin line) is the boundary condition at x = xq for (j3.86p . 
in which Og = f3g = 5g = Eg = 0, and 7^ = 1. The growth term G is modeled using the 
assumption that the resonant instability operates, i.e., G = Tj-esW [see Equation (|3.63p ]. 
where an artificial CR pressure spectrum was imposed, described by the expression 

P^,{x,p) = O.Spong— e-^i'^P-i^Po)' exp (- — ^] . (3.160) 
Po \ xo p J 

(this pressure was simulated and binned into the momentum and spatial grids in order to 
emulate the actual run, where the pressure Per is calculated by the Monte Carlo particle 
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transport routine). The corresponding solution given by (|3.1U3|) is: 



W{{),k)=W{xo,k)+ / Q.bpoul—e 
Jxo dx I po 



(Inp-lnpo)^ 




x' po\ p dx' 
xo p J\k uq 



W{xo, k) + Q.bpQU^ 




1 — exp 



k 



) 



(3.161) 



where ko = cBq/cpq = (muo/po) r^^ , and po = mpC. The result of the numerical integration, 
shown with the solid thick line, coincides perfectly with the analytic solution (|3.161|) shown 
with the triangular markers. 

The cascading term, proportional to 5g, along with the viscous dissipation in the 
term proportional to tested in the following two runs. 

In Figure 13.171 I illustrate the third test - the solution of (j3.86p with Og = Pg = 
and jg = 6g = Eg = 1; the growth rate, G, was chosen similarly to the previous example, 
but with p = lO^rripC; the cascading rate, IT, was taken in the form (j3.79p : and I chose the 
viscous dissipation model described by (|3.74p with k^ w 1.1 • lO'^rgo, corresponding to a 
temperature Tq = 10^ K in (|3.75l) . The resulting turbulence spectrum is consistent with the 
predictions of the Kolmogorov theory. The energy-containing interval of wavenumbers is 
around ko = eBo/cpo ~ 3- 10^^ r~Q, where the turbulence amplification takes place (see the 
previous example for the amplified spectrum not modified by cascading). Then follows the 
inertial interval, where the energy is carried from small k to the greater k by cascading; the 
power law index of the spectrum matches very well the Kolmogorov's k~^^'^ law described 
by the analytic solution (j3.117p . Finally, at short wavelengths, the dissipative interval is 
marked by the spectrum turning down exponentially due to the effect of viscous dissipation, 
L. It happens at /c ~ 0.1 k^. 

In another test of cascading, I confirm that, if the seed turbulence has a power-law 
form, and is not amplified, the cascading leads to the formation of an inertial interval with 
W oc k~^^^ followed by the dissipative interval, where the spectrum turns down exponen- 
tially. The setup of the run shown in Figure [3. 181 is similar to that of the previous example, 
but 7g = 0, and the seed turbulence spectrum contains more energy by a factor 10^ (the 
solid thin line) . The evolution with cascading leads to the formation of the spectrum shown 
with the thick solid line. Its slope is in agreement with the Kolmogorov's law indicated 
with the dashed line. 
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The tests presented above are only a few of the multitude of tests that I per- 
formed in order to confirm that my major contribution to the model, the magnetic field 
amplification module, adequately solves equation (j3.86p and calculates the effects of tur- 
bulence generation and dissipation on the flow. These effects are: plasma heating due to 
the turbulence dissipation [see equation (j3.76l and the text explaining it] , the contribution 
of turbulence to the momentum and energy balance, which affects the plasma flow (see 
Section r3.2.9p . and the determination of particle transport by the spectrum W{x,k) (Sec- 
tion 13. Sp . I should note that in several publications we used a model for magnetic field 
amplification that included the generation of waves traveling in both directions, but did not 
include cascading. This model and the corresponding numerical integrator are described in 
Appendix A. 

3.2.9 Turbulence and equations of motion 

The fundamental idea on which the Monte Carlo model, as well as simpler analytic 
models, is based, is that the dynamics of matter, particles and magnetic fields are described 
on scales much larger than the scale of turbulent fluctuations. That is, the model does not 
contain and describe the information about the spatial structure of stochastic flows and 
magnetic fields, substituting an averaged statistical description. This is expressed in the 
following approximations: 

• Instead of a field of turbulent fiuctuations of the plasma velocities, the model has the 
averaged flow speed u{x); 

• Instead of the spatial structure of magnetic fields B(r, t), the Fourier spectrum of 
fluctuations, averaged over a large enough volume surrounding a coordinate x, is 
used, denoted as W{x,k); 

• Instead of describing particle transport using the equations of motion based on the 
Lorentz force, the model employs a diffusion model, in which the mean free paths 
depend on W{x,k). This diffusion approach applies on scales on which the particles 
'lose memory' of their initial direction of motion, and these scales must be greater 
than the size of the turbulent structures scattering the particles. 
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Figure 3.15: Effect of flow compression on turbulence spectrum. 
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Figure 3.16: Amplification of turbulence spectrum. 
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Figure 3.17: Amplification and cascading of turbulence. 




Figure 3.18: Cascading of seed power law spectrum of turbulence. 
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The above approximations mean that the equations of motion describing u{x) 
must contain the properly averaged contributions of the turbulence to the fluxes of mass, 
momentum and energy. In this section, we present and explain these contributions. The 
equations and reasoning shown here are pertinent to the discussions in Sections 13.1.71 and 

One may calculate the flux of momentum and energy, accounting for the turbulent 
contribution, using the general expression for the energy density Wt, the stress tensor Tj^, 
and the energy flux q (e.g., equations (2.48), (2.49) and (2.67) in [120J) 

Wt = p(lu^ + -)+^, (3.162) 

Tik = PSik + pUiUk + —Sik — (3.163) 

OTT 47r 

/I n e P\ B X (u X B) 

Here 6ik is the Kronecker delta-symbol, and the index 't' in Wt indicates that this is the 
total energy density of the bulk flow, accelerated particles, and turbulence. 

For simplicity (also see the comment at the end of this section), let us assume that 
the spectrum of turbulence, W{x, k), is a power spectrum of Alfven waves traveling along 
the magnetic field Bq in a plasma moving at a constant speed uq with mass density po- 
Such waves induce perturbations of the matter velocity and magnetic field, and the total 
flow velocity u and total magnetic field B can be written as: 

U = Uq + 5Umexp[ik{x — {Ux ^ VA))i^t] = Uq + 5u (3.165) 

B = Bo + 5Bm exp [ik{x - (u^ =F VA))uJt] = Bq + 5B, (3.166) 

where (5B and 6u are the time and coordinate-dependent values of the fluctuations of the 
magnetic field and the plasma velocity in the wave, (5Bm and 5um are the amplitudes of 
these fluctuations, and Ux is the average x-component of the flow velocity, also denoted 
throughout this work as u. The it signs correspond to different polarization (i.e., directions 
of motion). For Alfven waves, the following properties must be listed: 5u = ziidB/y/AirpQ, 
Bq II Uq, 5B _L Bq, (5u _L Uq. Also, becausc Alfven waves are an incompressible motion of 
plasma, one may add these conditions: p = pQ, e = eo and P = Pq (here e is the internal 
energy, and P - the pressure of the gas). 
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Substituting the expressions for u and B from (|3.165|) and (|3.166|) into (|3.162|) . 
(|3.163p and (|3.164p . one may derive the quantities that interest us in the 1-D simulation: 
W{x,k), Txx and Qx- Note that these quantities have also been denoted above as <I>p and 

= po(huo + 5uf + ^) + ^^^^±^, (3.167) 



2 PoJ Stt 

'^P = Txx = Po + pul+ ^ ^ (3.168) 

^E = qx = pouo (l-{uo + 6uf + — + — ] + 
V2 Po PoJ 

^ (Bo + ^B)x[(uo + Hx(Bo + ^B)] 

Simplifying the vector operations and averaging over many wavelengths in x and many 
cycles in t (this leads to ((^B^) = 5B^/2 and {Su"^) = du'^/2), one gets: 

(Wt) = ^pong + eo + ^ + Qm^^ + ^) A (3.170) 

{<^p)^{Txx) = p«g + Po-^+(^)/2, (3.171) 
1 

{'^e) = (qx) = ^PoUq + (Po + eo)uo + 

+ (ipou„fai+ "°^^'";/°''< )/2. (3.172) 

In the following, we omit the averaging signs ( ). Associating the last terms in the above 
equations with the contributions of turbulence we have: 

W = Qm«^ + ^)a (3.173) 

= (^)a (3.174) 

-Pquq5u^^ + — 1 /2. (3.175) 



103 



Now, using the 'equipartition' characteristic of Alfven waves, i.e., the identity du^ = 
(5i?^/(47r/)o), and the definition of Alfven velocity va = -Bq/V^'^PO) we arrive at: 

W = W, + W^ = Y-^ + \^, (3.176) 
= \w, (3.177) 
= ^{uo^va)W. (3.178) 

In these equations, Wk = W„i are the energy densities of, respectively, kinetic and magnetic 
turbulent fiuctuations. Equations (|3.177p and (|3.178p define the 'pressure' (i.e., flux of the 
x-component of momentum in the x-direction) and the energy flux (in the x-direction) of 
turbulence. These quantities should be added to the corresponding fluxes of particles in 
order to account for turbulence in the momentum and energy balance; in other words, in 
order to account for turbulence in the equations of averaged motion. 

Let us discuss the equations (|3.177p and (j3.178p deflned above. First of all, they 
only strictly apply to Alfven waves (but, thankfully, of arbitrary amplitude). Nonlinear in- 
teractions between high amplitude waves and particles may, as explained in earlier sections, 
lead to the turbulent behavior characterized by cascading and by significant changes in the 
geometry of magnetic fields and random plasma velocities, invalidating (|3.177p and (|3.178p . 
Also, even without the transition to turbulence, these equations do not rigorously apply to 
any waves other than Alfven. For instance, the short-wavelength harmonics generated by 
Bell's instability are not Alfvenic; one may show that for waves at A: = kc/2 (the peak of the 
growth rate), the balance between the kinetic and magnetic energy density of these waves, 
Wk and Wm, is Wm = 2>Wk-, as opposed to Wk = Wm for Alfven waves. 

In the absence of a more detailed model of turbulence evolution that describes the 
geometry and dynamics of stochastic motions and fields in the plasma, one cannot expect 
to significantly improve the calculation of and F^,. However, I argue that, as shown 
by the example of Bell's harmonics, different geometry or dynamics of turbulence may just 
lead to changes in the factors such as 1/2 and 3/2 in equations (j3.177p and (|3.178p . One 
may hope that this would be a minor change, where by 'minor' I mean a change by a 
factor of a few. This is as much certainty as one may expect to achieve without describing 
the spatial structure of turbulence with a PIC or MHD model. That approach, as we 
saw earlier, is extremely computationally expensive, especially for nonlinear shocks that 
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require a large spatial and temporal dynamic range, and I accept the equations (|3.177|) and 
(|3.178p in the model for the sake of achieving the designated goal of this work: studying the 
nonlinear structure of shocks undergoing efficient particle acceleration and strong magnetic 
field amplification. 
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3.3 Particle transport 

The problem of diffusive transport of charged particles in magnetized plasmas is 
fundamental for plasma physics. In collisionless plasmas typically found in astrophysics, 
this transport is generally turbulent diffusion as particles propagate in stochastic magnetic 
fields and the associated stochastic plasma motions. The question usually asked is, given 
the spectrum (or a more complete description - correlation tensors) of turbulence, find the 
diffusion coefficient of a particle with a certain momentum p. In this work I used several 
approximations of diffusion coefficients, as described below. Each of these approximations 
has a its own domain of applicability. 

3.3.1 Bohm diffusion limit 

Bohm diffusion was first observed for electrons in a magnetized laboratory plasma 
|76j , but the Bohm diffusion model is often applied in astrophysics due to its simplicity. The 
principal assumption is that the plasma is magnetized and turbulent, so that a particle's 
mean free path between strong deflections is equal to its gyroradius, 

Asohm = (3.179) 

Here p is the momentum of the particle, and B is the magnetic field in the plasma. The 
corresponding diffusion coefficient, assuming isotropic diffusion, is 

Z^Bohm = (3.180) 

where v is the speed of the particle corresponding to momentum p. Note that for non- 
relativistic particles (p = rav <C mc), -DBohm oc p^, and for ultra-relativistic ones (p ^ mc, 
t; = c), the scaling is -Dsohm oc y. 

The Bohm approximation is clear and intuitive. It features two most important 
dependencies: the diffusion coefficient increases with the particle momentum, p, and de- 
creases with the magnetic field B. This diffusion model rests on the assumption that B is 
rather strong: it confines the particle gyromotion to scales on which the field itself varies 
significantly (so that the diffusive character of motion is effectuated). 
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3.3.2 Resonant scattering by Alfven waves 

When a uniform field, Bq, exists in a plasma on scales much larger than the sizes of 
particle gyroradii and turbulent harmonics, and a train of low amplitude AB <^ Bq Alfven 
waves travels along this field, the mean free path of an energetic particle along the uniform 
field can be estimated as 

Ares = (3.181) 

IT J- 



where 
and 



(3.182) 



feres = (3.183) 

(see [126j or [82] ) . In expression (|3.18ip , the numerator of the second fraction is the gyrora- 
dius of the particle (pj_ is the component of the particle's momentum transverse to the field 
Bo), and the denominator T is, within a factor, the energy density of Alfven waves (per 
unit logarithmic waveband din A; = 1) normalized to the energy density of the underlying 
uniform field. The energy density W(}t) in (j3.182p is taken at the resonant wavenumber /cres 
defined by (j3.183p . When T approaches 1, the mean free path shrinks down to the particle 
gyroradius, and the Bohm limit is realized. Increasing T further takes this theory beyond 
its applicability limits. 

3.3.3 Diffusion in short scale turbulent fluctuations 

If the bulk of the turbulence energy is in small-scale harmonics with respect to the 
particle mean free path, then the motion of the particle is nearly ballistic, with frequent 
and small deflections from the stochastic Lorentz force. The collision length of such motion 
can be expressed ([38], see also [117] and [7l]) as: 

'\ 

(3.184) 



4 ]?(? 

Xss{x,p) = -^2~ 



vr 



oo 

^ f W{x,k) 




This corresponds to a mean free path in the small-scale field, Agg, given by the expression 
Agg = rgg//cor) where rgg = cpjeB^^ is the gyroradius of the particle with momentum -p in the 
effective small-scale field, -Bgg, and /cor is equal to the correlation length of the small-scale 
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magnetic field (see below for exact definitions). This relationship is easy to understand. 
Consider a thought experiment: an energetic particle with momentum p is propagating 
through a medium consisting of regions of scale /cod each of which contains a magnetic field 
with magnitude Sgs, pointing in a different random direction in each region. In the course 
of the path Ass ^ ^cor; the particle encounters N = Ass/^cor ^ 1 such regions, and in each 
of them its momentum gets a random scattering in the amount Apv ^ F/S.t = eBsJcor/c 
(here F is the magnitude of Lorentz force, and At is the time of the particle crossing the 
region). Considering this process a random walk in p-space, the mean square deflection 
of momentum along the path Ass is (Ap)^ = A^(Apv)^ = Ks/lcor{eBsslcor/c)^ , and setting 
(Ap)^ = p^, corresponding to Ass being the mean free path, one can solve this equation 
to find Ass = (cp/eBss)'^ /Icor = ^"ss/^cor- This mean free path depends on the particle 
momentum as Ass oc as opposed to the Bohm behavior oc p, which is a significant 
difference. 

3.3.4 Low energy particle trapping by turbulent vortices 

Suppose the turbulence has a power law spectrum that contains a significant frac- 
tion of energy in the smallest scales (such a spectrum may be produced by cascading as 
described in Section [3. 2. 4p . A particle with a low enough energy will be effectively confined 
by resonant scattering on the small scale turbulence fluctuations. But its transport on 
scales greater than the correlation length of the turbulence (i.e., greater than the largest 
turbulent harmonics), which is of interest for the Monte Carlo code, may be significantly 
different from the directly applied model of resonant scattering transport. The efficient 
resonant scattering effectively confines the particles to the large-scale turbulent structures, 
and their diffusion on large scales is determined by the motions of the turbulence rather 
than the particles' own motion. A theoretical description of such transport is described by 
Bykov and Toptygin in [28], [29] and |120j . A rough approximation of their result is that, if 
the mean free path of a low energy particle due to resonant scattering is Arcs ^ ^cor , where 
/cor is the correlation length of the turbulence, then the diffusion coefficient of such particle 
on scales greater than /cor is on the order of 

D « Ujcor, (3.185) 
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where Uc is the typical speed of turbulent motions with correlation length /cor- This applies 
when D ^ vX^es, where v is the speed of the particle, and Ares is its mean free path between 
the resonant scatterings, meaning that the 'convective' diffusion coefficient (|3.185p is much 
greater than the resonant scattering coefficient. This situation is analogous to the convective 
diffusion of cream in a coffee cup. Pour the cream into the coffee and, even without stirring, 
it will spread through the cup in minutes. If one naively assumes molecular diffusion and 
estimates the time it takes the cream to diffuse from one end of the cup to another, this 
time will be much longer, on the order of hours. The discrepancy is successfully explained 
with a model similar to (I3.185P : molecules of the admixture are confined to the turbulent 
vortices in the medium (in the coffee cup, those are induced by the temperature difference 
between the top and the bottom, and by the energy introduced during the pouring of the 
coffee into the cup and of the cream into the coffee), and the propagation of the admixture 
is determined by the motion of these vortices rather than of the admixture with respect to 
the vortices. 

Another possibility of particle trapping in turbulent structures is when there is 
no short-scale turbulence to produce effective resonant scattering, but a particle has a low 
enough energy so that its gyroradius in the large scale turbulent magnetic field, r^, is small 
compared to the correlation length of the turbulence, /cor- Then the particle will gyrate 
around the turbulent magnetic fields, losing the memory of its initial direction of motion 
on the length comparable to /cor- If the particle's speed v ^ Uc (so that the turbulence is 
essentially stationary for the particle), then one may estimate the coefficient of diffusion of 
the particle on scales greater than /cor as 

D ^ v/cor, (3-186) 

or the effective mean free path of the particle as 

A « /cor- (3-187) 

This means that particles trapped in the turbulent vortices by gyration in the turbulent 
large-scale magnetic fields have a mean free path nearly independent of the particle energy 
and equal to the size of the turbulent vortices. More realistic models of this process may 
be necessary, because effects such as drifts in magnetic fields and time dependence of the 
vortex structure may change the dependence of the mean free path on the particle energy. 
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The above approximation applies to particle transport on scales greater than the 
turbulence correlation length. The transport of very low energy particles on smaller scales 
depends on geometry and evolution of the turbulent structures, which is beyond the reach 
of our model. The motion of magnetic field lines (sometimes called magnetic field line 
wandering) may be non-diffusive on small spatial scales, resulting in CR transport that 
cannot be described as diffusion (e.g., |101j ). 

3.3.5 Implementation of diffusion models in the Monte Carlo code 

Based on the theoretical models of particle transport outlined above, I imple- 
mented the corresponding mean free path prescriptions into the Monte Carlo code. When 
the model is run, the user can specify which prescription is to be used in the simulation. 
It allows the application of transport models of various degrees of physical accuracy and 
applicability to study their effects on the self-consistent shock structure. 

Bohm diffusion 

If the user specifies the Bohm regime of diffusion in the simulation, then given the 
momentum of the particle, p, measured in the plasma frame, the code will calculate the 
mean free path ABohm using (j3.179p . where for B it substitutes the effective local magnetic 
field, i?eff) defined in (j3.6ip . 

This is the simplest method of describing diffusion in the presence of efficient 
MFA. It should give an accurate (within an order of magnitude) estimate of the collision 
mean free paths for moderate energy cosmic rays. For the highest energy cosmic rays, 
with gyroradii greater than the magnetic field correlation length, the turbulence acts as 
small-scale magnetic fiuctuations, an Bohm diffusion is an overestimate of the confinement 
strength. For the lowest energy CRs and thermal particles, Bohm diffusion is also not a 
good approximation, because the particles may be trapped in magnetic structures, in which 
case their diffusion is determined by the evolution of the small scale turbulence rather than 
their own motion. 
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Resonant scattering 

I have the option of describing the particle scattering with a form similar to (j3.18ip 
in the simulation. If this model is adopted, it calculates the resonant wavenumber as given 
by (|3.183p . except that it uses the total momentum p instead of , and to calculate the mean 
free path, it uses (I3.18ip . but with p instead of p±. This replacement of the components 
of the particle momentum with its magnitude is done in order to account for the strong 
nature of the turbulence. Actually, if AB ^ Bq, then ()3.18ip is not applicable in all 
rigorousness, but I use this theory in order to grasp the most important qualitative behavior 
of the turbulent transport: the stronger the turbulent structures of scales comparable to 
the particle gyroradius, the more efficient is particle scattering. 

If the turbulence spectrum has the shape W = Wo{k / ko)^^ , which will hereafter 
be called the Bohm spectrum, then (j3.18ip gives a mean free path similar to the Bohm 
prescription (|3.179p . Namely, when ^/AirWokQ = Bq, and = 1, the two models match 
within a factor of 4/7r. The latter condition is equivalent to the condition that a unit 
logarithmic waveband dink = 1 contains the same amount of turbulent energy as the 
underlying magnetic field Bq. 

Thus, for relativistic particles. Ares oc p in a Bohm spectrum W{k) oc k^^. Steeper 
spectra of turbulence {W oc A:~^ for g > 1) give weaker dependencies of Ares on p. The 
spectrum W{k) oc k~'^ gives a constant Ares(p)- 

Hybrid model of diflfusion in strong turbulence 

It is useful to re-write equation (j3.184p as 

Ass(x,p) = f^, (3.188) 

^cor 

where is the particle gyroradius in the effective magnetic fields of the short-scale magnetic 
perturbations. In the model, I adopt the prescription (j3.188p . and generalize it with two 
assumptions, as outlined below, so it can be applied to particles of lower energies as well. 
The first assumption is that for a particle of momentum p, the local turbulence spectrum 
can be divided into the large-scale and the short-scale part, the wavenumber A;* being the 
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boundary between them. The effective large-scale magnetic field is then 



^i^l^ = ^ + \f Wi^^ k') dk\ (3.189) 
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the effective small-scale field is 

Bl{x,k,) _ 1 
Svr 2 

k 

and the correlation length of short-scale field /cor can be estimated as 



j W{x,k') dk\ (3.190) 



f?°W{x,k')/k' dk' 
/~W^(x,fc')d)t' • ^ ^ 

I define k^ using the condition rg[B\s)k^ = 1, where rg{B\s) = cp/eB\^ is the gyroradius of 
the particle in the large-scale magnetic field B\s. The latter is dependent on /c*, therefore 
a nonlinear equation must be solved at every point in space for every particle momentum 
in order to determine A;*. The second assumption is that the calculated \ss{p) does not 
increase as momentum p decreases. 

Let us comment on the physics behind the assumptions outlined above. Equa- 
tion ()3.188p applies if the turbulence is predominantly short scale. However, if the turbu- 
lence spectrum incorporates a wide range of wavenumbers (for example, the assumed up- 
stream spectrum W oc fc^^) then a good quasi-linear approximation to the particle transport 
properties is the resonant scattering prescription (j3.181|) (see, e.g., [51 1122) and references 
therein). However, for a spectrum similar to (j3.60p and k^ <^ k^i^^, the mean free path 
(|3.188p can be represented after some simple mathematical transformations as 

eB,A^kwlx,ky ^^-^^^^ 
where k = eB\s/{cp). For B\s w Bq (weak turbulence case), this is precisely the resonant 
scattering mean free path (|3.18ip . and for Bis > (strong perturbations), it may be a 
good generalization of the latter. Therefore, dividing the turbulence spectrum at /c* allows 
one to correctly describe the mean free path of intermediate-energy particles using (j3.188p . 
along with the high energy particles. 

The second assumption, that of monotonic behavior of A(x,p) with respect to 
p, doesn't influence the case of a power-law turbulence spectrum, but affects the diffusive 
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transport of low energy particles in case of turbulence with a marked concentration of energy 
around a wavenumber k^, i.e. containing strong vortices of size \/k^. Indeed, assume for 
simplicity a Gaussian spectrum W{x,k) oc exp [(/c — fcv)^/(2o"^)] , where a is the width of 
the spectrum. If the particle momentum p is large enough so that k^, < fc^, and r^^ ^ /cor 
holds, then the particle is scattered by frequent deflections in the short-scale magnetic field 
of the vortices, and equation (j3.188p applies unconditionally. However, a particle with a 
low enough momentum so that k^ ^ k^ will find itself trapped in the large-scale magnetic 
fields of the vortices, and one may assume that its transport on scales larger that 1/k^ is 
diffusive (see equation (j3.187p . with the effective mean free path 

\K,l/k^. (3.193) 

Now consider the above Gaussian spectrum. The prescription (|3.188p . with r^s and /cor 
determined using the k^ formalism, does not describe the trapping of the low energy parti- 
cles. However, at such momentum ptr that k^ ~ k^ for this momentum, the magnetic field 
i?is ~ ^ Bq, (assuming strong turbulence), and lowering the value of p will lead to an 
exponentially rapid decrease of B^s, and an equally rapid increase in rgg, which will make 
A = rgg/Zcor unphysically increase for smaller p. The monotonicity assumption will correct 
this unphysical behavior by fixing \{x,p) at the value A(x,ptr) for p < ptr- And this value 
will approximately be (j3.193p . because ptr corresponds to ~ 1/A;v ~ /cor- 

Summarizing, I state that I choose the mean free path of particles with momentum 
p according to (|3.188p . where rgs and /cor are calculated for the short-scale part of the mag- 
netic field, k > k^, and force this prescription to be monotonic in p for low momenta. The 
reasoning provided above shows that our prescription properly describes particle transport 
a) for high p particles in short-scale field, as per the derivation of (j3.188p : b) for intermediate 
to low p in a power-law turbulence spectrum, assuming resonant scattering, and c) for low 
energy particles in large-scale turbulent vortices, assuming particle trapping. In between 
these important regimes, the prescription provides an interpolation. 
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3.4 Parallel computing with MPI 

The code used for this dissertation is written for parallel processing using the MPI 
(Message Passing Interface) protocol. In this section I will summarize the parallelization 
algorithm, outline its advantages and drawbacks, and present a performance test. 

By far the most time consuming part of the simulation is the Monte Carlo transport 
of particles that simulates the Fermi-I acceleration process. This procedure is intrinsically 
very well suited for parallel computing, because particles are propagated one after another, 
and each particle's motion within an iteration is completely independent of any other par- 
ticle's history. Multiple particles are only required in order to decrease random deviations 
of the results, i.e. to 'improve statistics'. This also means that the quality of random 
numbers is not a major issue of concern for this Monte Carlo code, because even if the 
random numbers are correlated within a sequence, correlated between different processors, 
or do not continuously fill their range of definition, it does not affect the quality of results. 
That is because the trajectory of each particle depends not only on the latest scattering 
outcome, but also on the previous history of acceleration of this particle, which effectively 
diminishes any possible correlations in particle histories due to the imperfections of the 
random numbers used by the Monte Carlo codq^. 

I implemented the following algorithm of parallelization of the calculations. First, 
a 'master' processor divides the user-specified number of particles equally between the 
available processors, including itself. Then each processor (the 'slaves' and the 'master') 
performs one iteration, i.e. propagates the particles it is responsible for until they reach the 
highest achievable energies, and the iteration terminates. After a processor completes its 
iteration, it returns the output, (the particle distribution function /(p) and its moments) to 
the 'master' processor (which also performs its iteration equally with the other processors, 
and returns the collected information to itself). The 'master' processor then averages the 
incoming results (which improves the statistical certainty of the calculated particle distribu- 
tion, of momentum and energy fluxes, and of the increments of field-amplifying instabilities) 
and uses them to calculate magnetic field amplification and precursor smoothing. These 
procedures are not easily parallelized, but they take relatively little time, and I chose to 



®The random number generator used in the code is an excellent match for the single-processor version of 
the Monte Carlo simulation, but was not specifically designed for parallel processing. This, however, turns 
out not to be a problem for the reasons stated in the text. 
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leave them to just one processor. After that, the 'master' processor gives the other proces- 
sors the updated flow speed u{x), the re-iterated magnetic turbulence W{x,k), and each 
processor computes the corresponding mean free path prescription X{x,p) and performs 
another iteration. This cycle continues until the self-consistent solution is derived. 

The primary advantage of this procedure is its ultimate simplicity. In terms ac- 
cepted in the parallel computing field, this is an 'embarrassingly parallel' code, which means 
that the interactions between processors take place very infrequently (in practice, they ex- 
change several megabytes of data once every several minutes). Another advantage is that 
one processor's runtime performance does not affect another processor's particle history. 
It is a welcome feature of the method, because it makes it easier to debug, if problems 
arise: the results, including the run-time errors, are reproducible. We must note here that 
the sequences of random numbers generated by the code are, in fact, deterministic: in two 
identical runs executed at different times, the random number sequences and the final re- 
sults will be identical. The same applies to the version of the code with parallel computing 
performed as described above. 

A disadvantage of this method is that there may be situations when most proces- 
sors had finished their iterations, but must wait for one processor working on a particle with 
an 'unfortunate', long history of acceleration. Computing time is lost in this case, because 
the duration of every Monte Carlo iteration is as long as the worst processor's performance. 
This is not a major issue of concern when the number of particles per processor is large, 
but when many processors are available, and each gets only a few particles, the deviation 
of the worst processor's performance from the average performance may be significant. 

In Table [312] I listed the results of a simulation similar to that done in Section TS.l.Sl 
I executed 5 runs, with identical input parameters, but with different numbers of processors: 
1, 2, 4, 8 and 16. Each run obtained a self-consistent shock structure, and the results were 
identical in all runs within small statistical deviations. 

Column 'A^proc' lists the number of processors used in the run. There were a total 
of Np = 160 particles in each ruj^, and they were equally divided between the processors, 

^This is not the number of thermal particles. A numerical procedure called 'particle splitting' is used 
in the Monte Carlo model, which allows to maintain nearly equal number of particles at any energy - a 
necessary condition to simulate rapidly decreasing particle spectra over many decades of the energy. I did 
not describe the 'particle splitting' in this dissertation, because it is purely technical, and was explained in 
the literature (e.g., [72]). 
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Table 3.2: Test of performance boost with parallel computing 



N 

^ ' proc 


-^p/-^proc 


Time, s 


Speedup 


1 


160 


7622 


1.0 


2 


80 


5614 


1.4 


4 


40 


3049 


2.5 


8 


20 


2092 


3.6 


16 


10 


1688 


4.5 



as shown in column 'Np/Nproc. The time that the each simulation took is listed in seconds, 
and the column 'Speedup' shows the ratio of the time of the iteration with 1 processor (i.e., 
without parallel computing) to that of the parallelized run. 

The results show that the parallelization does lead to an increase in the speed of 
the calculations, but the speedup is proportional to, approximately, {Nproc)'^'^, which is not 
very efficient. An alternative parallelization algorithm that may improve the efficiency is 
reserving the 'master' processor for dispatching particles between the 'slaves' in the run- 
time, i.e., each 'slave' gets a new particle from the 'master' as soon as it finishes with the 
previous one. This way the situation when many processors await a few 'unlucky' ones to 
finish with their iterations will not last as long. 

However, because this procedure is incompatible with the reproducibility of results 
(unless each particle has its own random number sequence, that is stored and passed between 
the processors; care must be taken in this case to ensure that correlations between particle 
histories do not occur when particle splitting is performed). This makes it difficult to debug 
the code (see above), I chose to stay with the currently implemented, less efficient, but more 
predictable scheme. Despite the less than perfect scaling of performance with the number 
of processors, it allows me to achieve reasonable computation times even with the available 
modest computational resources (8-16 processors per run). Typical run times are seconds 
to minutes for test runs, and around one day for a self-consistent simulation with a realistic 
dynamic range and small enough statistical deviations. 
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Chapter 4 

Applications of the model 

In this chapter I will show the basic results of the model of magnetic field ampli- 
fication in collisionless shocks based on the Monte Carlo simulation of DSA. Some of these 
results have appeared in peer-reviewed publications (Sections 14. H 14.2] and 14.3]) . and some 
will soon be submitted for publication (section 14. 4p . 

For the published material, in this Chapter I only provide a condensed version of 
the articles. For the work that has not yet appeared in press (Sections 14.51 t^is reader 
will find an outline of the proposed direction of research, a presentation of the preliminary 
results and a discussion of the applicability to astrophysical research. 
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4.1 Turbulence growth rate and self-consistent solutions 

In 1122]^ . we introduced a Monte Carlo model of nonlinear diffusive shock accel- 
eration allowing for the generation of large-amplitude magnetic turbulence, i.e., AB S> Bq, 
where Bq is the ambient magnetic field. The model is the first to include strong wave gen- 
eration, efficient particle acceleration to relativistic energies in nonrelativistic shocks, and 
thermal particle injection in an internally self-consistent manner. In order to describe the 
field growth rate in the regime of strong fluctuations, we use a parameterization that is con- 
sistent with the resonant quasi-linear growth rate in the weak turbulence limit. We believe 
our parameterization spans the range between maximum and minimum rates of fluctuation 
growth. 

We find that the upstream magnetic field Bq can be amplified by large factors and 
show that this amplification depends strongly on the ambient Alfven Mach number. We 
also show that in the nonlinear model large increases in B do not necessarily translate into 
a large increase in the maximum particle momentum a particular shock can produce. The 
most direct application of our results will be to estimate magnetic fields amplified by strong 
cosmic-ray modified shocks in supernova remnants. 

4.1.1 Model 

In |122] . we described the amplification of magnetic turbulence by the following 
set of equations: 

U+ + U^ ^ dx 



+ — {U.- U+) , (4.2) 

which were solved iteratively in the MC simulation. This system describes the development 
of the resonant cosmic ray streaming instability of Alfven waves along with the processes 
^The results presented here first appeared in [122] and largely are reproduced from this publication. 
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of wave amplitude increase due to the plasma compression, and of interactions between 
waves traveling in opposite directions (see Section [3.2.41 and Appendix A). For the growth 
of Alfven waves in quasi-linear theory, Vg = va, where va = Bo/y^Airp^x) is the Alfven 
speed calculated with the non-amplified field and p{x) is the matter density at position x. 
This choice of Vg provides a lower limit on the amplification rate for the nonlinear regime, 
AB ^ Bq, and was used in |3J. If, on the contrary, we define Vg using the amplified field, 
i.e., Vg = Be{i{x) / \J Airp^x) , it reflects the situation where the growth rate is determined by 
the maximum gradient of Pct{x,p) along the fluctuating field lines. This provides an upper 
limit on the wave growth rate and was used in [llj. The real situation should lie between 
the two extremes for Vg- For this preliminary work, we vary Vg between the two limits, 
i.e., we introduce a parameter, < /aif < 1, such that 

/alf} , (4.3) 

and Vg varies linearly between va (for /aif = 0) and Bcs / \/ 4:iTp{x) (for /aif = 1). 

Finally, we assume a Bohm model for diffusion. The mean free path of a particle 
with momentum p at position x is taken to be equal to the gyroradius of this particle in the 
amplified field, i.e., \{x,p) = rg{x,p) = pc/[qBcs{x)], and the diffusion coefficient is then 
D{x,p) = Xv/3, where v is the particle speed, and B^s was defined according to (j3.6ip . 

The Monte Carlo code used here was the original simulation developed by Ellison 
and co-workers, not the version developed by the author of this dissertation for the problem 
of magnetic field amplification. I have confirmed that the latter model reproduces the results 
presented here. 

4.1.2 Results 

In all of the following examples we set the shock speed uq = 5000 km s~^, the 
unshocked proton number density UpQ = lcm~^, and the unshocked proton temperature 
To = 10^ K. For simplicity, the electron temperature is set to zero and the electron contri- 
bution to the jump conditions is ignored. With these parameters, the sonic Mach number 
Ms ~ 43 and the Alfven Mach number Maif ~ 2300(l/xG/5o). 



Vg = va<1 + 



Bn 
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With and Without Magnetic Field Amplification 

Figure [d] shows self-consistent solutions for four shocks, obtained with /aif = 0. 
Note that the horizontal scale has units of rg{uo) = mpUo/(ei?o) and is divided at x = 
—5rg{uo) between a linear and log scale. The the heavy dotted curves show results without 
amplification and all other curves are with amplification. The heavy solid and dotted curves 
have dpEB = — lO^rg(no), the dashed curve has dpEB = — lOOOrg(Mo)) the light solid curve 
has dpEB = -lO^rg(uo). 

First, we compare the results shown with heavy-weight solid curves to those shown 
with heavy-weight dotted curves. The heavy solid curves were determined with S-field 
amplification while the dotted curves were determined with a constant B(,f[{x) = Bq. All 
other input parameters were the same for these two models, and an upstream free escape 
boundary was placed at dpEB = — lO^rg(uo), where rg{uo) = nipUoc/ {eBo). The most 
striking aspect of this comparison is the increase in i?cfT(a;) when field amplification is 
included (bottom panels). The magnetic field goes from Bcs{x — oo) = 30 fiG, to B^s > 
1000 /uG for X > 0, and this factor of > 30 increase in B will influence the shock structure 
and the particle distributions in important ways. The solution without B-field ampliflcation 
(dotted curves) has a considerably larger rtot than the one with amplification, i.e., for no in- 
field amplification, rtot — 22, and with i?-field amplification (heavy solid curves), rtot — nEl 
This difference in overall compression results because the wave pressure Pw is much larger in 
the field amplified case making the plasma less compressible. A large rtot means that high 
energy particles with long diffusion lengths get accelerated very efficiently and, therefore, 
the fraction of particles injected must decrease accordingly to conserve energy. The shock 
structure adjusts so weakened injection (i.e., a small r^uh) just balances the more efficient 
acceleration produced by a large rtot- Since r^nh largely determines the plasma heating, the 
more efficiently a shock accelerates particles causing rtot to increase, the less efficiently the 
plasma is heated. 

In Figure we show the phase space distributions, f{p), for the shocks shown in 

Figure HTTJ These spectra are multiplied by [p/(?Tipc)]^ and are calculated downstream from 

the shock in the shock rest frame. For the two cases with the same parameters except field 

^See [12] for a discussion of how very large rtot's can result in high Mach number shocks if only adiabatic 
heating is included in the precursor. The uncertainty on the compression ratios for the examples in this 
paper is typically ±10%. 
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Figure 4.1: Shock structure with and without MFA 
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Figure 4.2: Phase space distributions with and without MFA 



122 



amplification, we note that the amphfied field case (heavy solid curve) obtains a higher Pmax 
and has a higher shocked temperature (indicated by the shift of the "thermal" peak and 
caused by the larger r^uh) than the case with no field amplification (heavy dotted curves). 
It is significant that the increase in Pmax is modest even though B increases by more than a 
factor or 30 with field amplification. We emphasize that Pmax as such is not a free parameter 
in this model; Pmax is determined self-consistently once the size of the shock system, i.e., 
c^FEBi and the other environmental parameters are set. 

In order to show the effect of changing dpEB; we include in Figs. 14.11 and 14.21 field 
amplification shocks with the same parameters except that dpEB is changed to — lOOOrg(no) 
(dashed curves) and —10^ rg{uQ) (light-weight solid curves). From Figure HT2l it's clear that 
Pmax scales approximately as dpEB and that the concave nature of f{p) is more pronounced 
for larger Pmax- The field amplification also increases with Pmax, but the increase between 
the dpEB = — lOOOrg(uo) and dpEB = —10^ rg{uo) cases is less than a factor of two (bottom 
panels of Figure 14. ip . 

In Figure Sis] we show the energy density in magnetic turbulence, C/+(x,A;) + 
f/_(x,fc), the diffusion coefficient, D{x,p), and particle distributions as functions of k and 
p at three different positions in the shock. All of these plots are for the example shown 
with dashed curves in Figs. 14.11 and 14.21 (i.e., with dpEB = — lOOOrg(uo)). The solid curve 
is calculated downstream from the shock, the dashed curve is calculated at x = —rg{uo) 
upstream from the subshock, and the dotted curve is calculated at a; = — lOOrg(tio) upstream 
from the subshock. 

The efficiency of the shock acceleration process can be inferred from Figure 14.41 
It shows the number density of particles with momentum greater than p, i.e., A^(> p), the 
energy density in particles with momentum greater than p, i.e., E{> p), for the shocks shown 
in Figs. 14.11 and 14.21 with heavy solid and dotted curves. The plots in Figure 14.41 indicate 
that the shocks are extremely efficient accelerators with > 50% of the energy density in f{p) 
placed in relativistic particles (i.e., p > mpc). The actual energy efficiencies are considerably 
higher since the escaping particles carry away a larger fraction of the total energy than is 
placed in magnetic turbulence. With Qcsc included, well over 50% of the total shock energy 
is placed in relativistic particles. Despite this high energy efficiency, the fraction of total 
particles that become relativistic is small, i.e., A^(> p = rupc) ~ 10~^ in both cases. 
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Figure 4.3: Turbulence and particle spectra with MFA 
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Figure 4.4: Acceleration efficiency with and without MFA 
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The effect of magnetic field amplification on the number of particles injected is 
evident in the left-hand curves. The larger rgub (solid curve) results in more downstream 
particles being injected into the Fermi mechanism with amplification than without. While 
it is hard to see from Figure I4.4[ when the escaping energy flux is included, the shock with 
-B-field amplification puts a considerably smaller fraction of energy in relativistic particles 
than the shock without amplification. Again, injection depends in a nonlinear fashion on 
the shock parameters and the subshock strength will adjust to ensure that just the right 
amount of injection occurs so that momentum and energy are conserved. 

Alfven Mach Number Dependence 

In Figure 14.51 we show three examples where Bq was varied and all other input 
parameters were kept constant. 

In all panels, the solid curves are for Bq = 0.3 /iG, the dashed curves are for 
Bq = SfiG, and the dotted curves are for Bq = SOfiG. The FEB is placed at the same 
physical distance in all cases with dpEB = — 1.7x10^^ m. Note that B^g increases most 
strongly for Bq = 0.3 /uG, but that the pressure in magnetic turbulence (bottom panel) 
does not gets above ~ 10% of the total pressure, which but still contains a significant 
fraction of the total pressure. 

The resulting overall compression ratios are: rtot — 9 for Bq = 0.3 ^G, rtot — 12 
for Bq = SfiG, rtot — 8 for Bq = 30 ^G, values consistent, within statistical errors, with 
Qcsci as indicated in the energy flux panels. As for the self-consistent amplified magnetic 
fields in these examples, B2/BQ ~ 400 for Bq = 0.3 /iG, B2/BQ ~ 150 for Bq = 3 fiG, and 
B2/BQ ~ 30 for Bq = 30 fiG. 

Wave Amplification factor, /aif 

All of the examples shown so far have used the minimum amplification factor 
/aif = (equation 14. 3p . We now investigate the effects of varying /aif between and 1 so 
that Vg varies between Va{x) and Bcs{x) j \J Attp(x). The other shock parameters are the 
same as used for the dashed curves in Figure SHI i.e., uq = 5000 km s^^, Bq = 30 /xG, and 
dpEB = -lOOOrg(iio). 

Figure shows u{x)/uq and Bf,Q{x) / Bq for /^if = 0, 0.1, 0.5, and 1 as indicated. 
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Figure 4.5: Comparison of shocks with different far upstream fields Bq. 
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Figure 4.6: Shocks with varying /aif 
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In all cases, mq = 5000km s^^, Bq = 30 ^G, and dpEB = — lOOOrg(no). The top panels show 
that increasing the growth rate (increasing /aif and therefore Vq) produces a large change 
in the shock structure and causes the overall shock compression ratio, rtot, to decrease. 
The decrease in rtot signifies a decrease in the acceleration efficiency and a decrease in the 
fraction of energy that escapes at the FEB, and the subshock compression adjusts to ensure 
conservation of momentum and energy. It is interesting to note that Tsub increases as rtot 
decreases and becomes greater than 4 for /aif > 0.5. In contrast to the strong modification 
oiu{x), there is little difference in Bcs[x) / Bq (bottom panels of Figure l4.6p and little change 

(Figure [471), between these examples. The fact that increasing the wave growth 
rate decreases the acceleration efficiency shows the nonlinear nature of the wave generation 
process. The most important reason for this is that the magnetic pressure P^, becomes 
significant compared to p{x)v?'{x) when /aif 1. The wave pressure causes the shock to 
be less compressive overall and forces rtot down. 

In Figure HTTl we show the distribution functions for the four examples of Figure IT6l 
and note that the low momentum peaks shift upward significantly with increasing /aif. As 
we have emphasized, the injection efficiency, i.e., the fraction of particles that enter the 
Fermi process, must adjust to conserve momentum and energy and the low momentum 
peaks shift as a result of this. The solid dots in Figure 14.71 roughly indicate the injection 
point separating "thermal" and superthermal particles for the two extreme cases of /aif = 
and 1. The first thing to note is that this injection point is not well defined, a consequence 
of the fact that the MC model doesn't distinguish between "thermal" and "nonthermal" 
particles. Once the shock has become smooth, the injection process is smooth and the 
superthermal population smoothly emerges from the quasi-thermal populationlf] Neverthe- 
less, the approximate momentum where the superthermal population develops, pinj, can be 
estimated and we mark this position with solid dots for /aif = and 1. What is illustrated 
by this is that the injection point shifts, relative to the post-shock distribution, when /aif is 
varied. This implies that, if injection is parameterized, the parameterization must somehow 
be connected to modifications in the shock structure. 

^We note that the smooth emergence of a superthermal tail has been seen in spacecraft observations of 
the quasi-parallel Earth bow shock (i.e., [51]) and at interplanetary traveling shocks (i.e., [8]). 
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4.1.3 Discussion 

We have introduced a model of diffusive shock acceleration which couples thermal 
particle injection, nonlinear shock structure, magnetic field amplification, and the self- 
consistent determination of the maximum particle momentum. This is a first step toward a 
more complete solution and, in this preliminary work, we make a number of approximations 
dealing mainly with the plasma physics of wave growth. Keeping in mind that our results are 
subject to the validity of our approximations, we reach a number of interesting conclusions. 

First, our calculations find that efficient shock acceleration can amplify ambient 
magnetic fields by large factors and are generally consistent with the large fields believed 
to exist at blast waves in young SNRs, although we have not attempted a detailed fit 
to SNR observations in this paper. More specifically, we find that the amplification, in 
terms of the downstream to far upstream field ratio B2/BQ, is a strong function of Alfven 
Mach number, with weak ambient fields being amplified more than strong ones. For the 
range of examples shown in Figure 14.51 B2/BQ ~ 30 for Maif ~ 80 and B2 / Bq ~ 400 for 
-^aif ~ 8000. Qualitatively, a strong correlation between amplification and Maif should not 
depend strongly on our approximations and may have important consequences. Considering 
that evidence for radio emission at reverse shocks in SNRs has been reported (see }66j , for 
example) and the strong amplification of low fields we see here, it may be possible for reverse 
shocks in young SNRs to accelerate electrons to relativistic energies and produce radio 
synchrotron emission. If similar effects occur in relativistic shocks, these large amplification 
factors will be critical for the internal shocks presumed to exist in 7-ray bursts (GRBs). 
Even if large i?-field amplification is confined to nonrelativistic shocks, amplification will 
be important for understanding GRB afterglows, in the stages when the expanding fireball 
has slowed down. 

As expected, amplifying the magnetic field leads to a greater maximum particle 
momentum, Pmax, a given shock can produce. Quantifying Pmax is one of the outstanding 
problems in shock physics because of the difficulty in obtaining parameters for typical 
SNRs that allow the production of cosmic rays to energies at and above the CR knee 
near 10^^ eV. Assuming that acceleration is truncated by the size of the shock system, we 
determine Pmax from a physical constraint: the relevant parameter is the distance to the 
free escape boundary in diffusion lengths. Our results show that Pmax does increase when 
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field amplification is included, but the increase is considerably less than the amplification 
factor at the shock B2/BQ (compare the heavy dotted and heavy solid curves in Figure [12]) • 
The main reason for this is that high momentum particles have long diffusion lengths, and 
the weak precursor magnetic field well upstream from the subshock determines Pmax- If 
the shock size, in our case dpEB; limits acceleration, Pmax will be considerably less than 
crude estimates using a spatially independent B2 (see also Section [42]) . On the other hand, 
particles spend a large fraction of their time downstream from the shock where the field 
is high and collision times are short. If shock age limits acceleration rather than size, we 
expect the increase in Pmax from the amplified field to be closer to the amplification factor, 
B2/BQ. 

Finally, it is well known that DSA is inherently efficient. Field amplification re- 
duces the fraction of shock ram kinetic energy that is placed in relativistic particles but, 
at least for the limited examples we show here, the overall acceleration process remains 
extremely efficient. Even with large increases in i?efr(a;), well over 50% of the shock en- 
ergy can go into relativistic particles (Figure [4. 4p . As in all self-consistent calculations, the 
injection efficiency must adjust to conserve momentum and energy. In comparing shocks 
with and without field amplification, we find that field amplification lowers rtot and, there- 
fore, individual energetic particles are, on average, accelerated less efficiently. In order to 
conserve momentum and energy, this means that more thermal particles must be injected 
when amplification occurs. The shock accomplishes this by establishing a strong subshock 
which not only injects a larger fraction of particles, but also more strongly heats the down- 
stream plasma. This establishes a nonlinear connection between the field amplification, the 
production of cosmic rays, and the X-ray emission from the shocked heated plasma. 
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4.2 Impact of MFA on the maximum particle energy 

Evidence is accumulatinj^ suggesting that collisionless shocks in supernova rem- 
nants (SNRs) can amphfy the interstellar magnetic field to hundreds of microgauss or even 
milligauss levels, as recently claimed for SNR RX J1713. 7-3946 j ll9j . Ironically, the evi- 
dence for large magnetic fields and, therefore, nonlinear MFA is obtained exclusively from 
radiation emitted by relativistic electrons, while the nonlinear processes responsible for 
MFA are driven by the efficient acceleration of relativistic ions, mainly protons. 

Here we address a single question: Can, as asserted by Uchiyama et al. [119J, the 
large amplified fields inferred for electrons from radiation losses in a nonlinear shock also 
determine the maximum proton energy produced in the SNR shock? We find the answer to 
be no because the inevitable nonlinear shock modification (due to efficient DSA) and the 
magnetic field variation in the shock precursor (due to MFA) make the maximum proton 
energy smaller than what is expected without accounting for these effects. 

Our result is similar to that found by [21j in a time-dependent calculation of DSA 
where the acceleration is limited by the age of the shock rather than the size, an indication 
that the nonlinear effects we discuss are robust. 

4.2.1 Model 

In a size limited shock, the proton maximum energy, E^^^, will be determined 
when the upstream diffusion length of the most energetic protons becomes comparable to 
the confinement size of the shock, typically some fraction of the shock radius. We model the 
confinement size with a free escape boundary (FEB) at a distance LpEB in front of the shock. 
Protons that reach this position stream freely away from the shock without producing any 
more magnetic turbulence. Therefore, for Bohm diffusion (see j53j for details), 

E^^^ oc LpEBnskSsk, (4.4) 

where is the upstream fiow speed. For a quasi-parallel, unmodified (UM) shock with 
no MFA, Ugki^sk = uqBq = uqB2, uq being the shock speed and B2 being the downstream 
magnetic field derivable from synchrotron emission of accelerated electrons. However, for 
a nonlinear (NL) CR modified shock of the same physical confinement size, LfeB; the 
*This section presents, in a condensed form, our publication [53) . 



133 



maximum proton energy -E^^^Inl will be determined by some mean value {u{x)B{x)) , 
giving 

(4.5) 



E^^^Inl _ {u{x)B{x)) 



^^^^luM U0B2 

For a strongly modified shock, {u{x)B{x)) <^ U0B2, and in the following we determine 
{u{x)B{x)) /{UQB2) using the Monte Carlo model described in detail in |122j . 

The Monte Carlo model we use (see [ 122] for full details) calculates NL DSA and 
the magnetic turbulence produced in a steady-state, plane-parallel shock precursor by the 
CR streaming instability. We self-consistently determine the nonlinear shock structure [i.e., 
u{x) vs. x], the MFA [Bq{({x) vs. x], and the thermal particle injection^ 

The NL results we investigate do not depend qualitatively on the particular shock 
parameters as long as the sonic Mach number is large enough to result in efficient DSA. 
Here, we use a shock speed = uq = 3000 km s^^, sonic Mach number Ms ~ 30, plasma 
density nisM = 1 protons cm~^, and Bq = B^sm = 10 /iG, yielding an Alfven Mach number 
Mj^ ~ 140. To these parameters we add a FEB boundary at Lfeb ~ 0.1 pc, correspond- 
ing to lO^Tgo, where rgo = mpUQc/{eBQ). This size is comparable to the hot spots in 
SNR RX J1713. 7-3946 and produces a proton energy ~ 10^^ eV in our unmodified shock 
approximation. Using the above parameters, we simulate two cases: a nonlinear solu- 
tion, where B is amplified from an upstream value Bq = 10 fiG to a downstream value 
B2 = 450 /iG (obtained self-consistently by our model), and an unmodified solution with a 
magnetic field set equal everywhere to B2 = 450 fiG. In these two cases we look at E^^^ to 
see how the prediction of the NL model, conserving momentum and energy, compares to 
the prediction of the UM model, implicitly assumed by |119j . The information about the 
maximum energy of electrons (which are not included in our calculations) can be inferred 
graphically from the plot of the acceleration time (see Fig. 14. 9p . 

4.2.2 Results 

Figure HiH] shows the shock structure, u{x), the effective magnetic field after am- 
plification, Bcsix), and u{x)Bcf[ix)/{uoB2), for the unmodified case (dashed lines), and the 
nonlinear case (solid lines). The bottom panel shows the energy flux, normalized to the far 



^Note that the Monte Carlo model ignores the dynamic effects of electrons and the NL shock structure is 
determined solely from the pressure of the accelerated protons and of the amplified magnetic fields. While 
electron acceleration can be modeled (e.g., [7]), we only show proton spectra here. 
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upstream value, for the NL case. The smoothing of u{x), the weak subshock (rsub — 2.9), 
and the increase in rtot above 4 (rtot — 9) are clearly present in the top panel for the NL 
case. These three effects must occur to conserve momentum and energy if CRs are efficiently 
accelerated. The quantity u{x)B(,s{x) / [uqB2) ~ 0.1 over most of the precursor in the NL 
case. 

In Figure [491 we show the momentum distributions functions, f{p) (multiplied by 
p^), and the acceleration time, race- The NL effects evident in Fig. 14.81 result in, 

Pj?r/« < 0.1 , (4.6) 

and a longer race to a given momentum. Here, p^f^ = ^^'^'^^Inl/c and ~ ^^'^^Ium/c 
We have not attempted a detailed fit to SNR RX J1713. 7-3946, but note that the concave 
shape of our proton spectrum, above the thermal peak, is similar to that obtained by [15] 
who find a good fit to the data, including the HESS TeV observations [68j. 

The UM result is obtained from the Monte Carlo simulation assuming the same 
"thermal leakage" model for injection as in the NL result (e.g., |72j). This injection scheme 
works self-consistently with modifications in the shock structure and overall compression 
ratio, rtot 5 to conserve momentum and energy in the NL case. In the UM case, the shock 
structure and rtot are not adjusted and the thermal leakage model produces far too many 
injected particles to conserve momentum and energy. For the UM shock to become a test- 
particle shock with energy conservation, far fewer particles would need to be injected so 
that the normalization of the superthermal f{p) oc power law becomes low enough, 
relative to the thermal peak, so that in contains an insignificant fraction of the total shock 
ram kinetic energy. Since we are only interested in comparing p™*^^ in the two cases, the 
normalization of the unmodified power law is unimportant since P\jyl only depends on Lfeb- 

4.2.3 Discussion 

The possibility of strong MFA in SNR shocks has been strengthened by the recent 
observations of rapid time variability in hot spots in SNR RX J1713. 7-3946 by |119j . If we 
accept the conclusions of [119], the ~ 1 yr variations in X-ray emission in some hot spots 
stem from radiation losses for electrons and indicate magnetic fields on the order of 1 mG. 
Such large fields would almost certainly be caused by MFA occurring simultaneously with 
the efficient production of CR ions in DSA. 
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Flow structure: unmodified versus nonlinear. 
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Figure 4.9: Proton spectra and acceleration times: unmodified versus nonlinear. 
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While a number of other interpretations of the X-ray and broadband emission 
in SNR RX J1713. 7-3946 have concluded that the magnetic field present in the particle 
acceleration site is considerably less than ImG (e.g., [52 | \8 ^ [T5 l llOOj ). we have shown that 
even if the magnetic field inferred from electron radiation losses is as high as |119] claim, 
the underlying physics of MFA in DSA shows that this field cannot be simply applied to 
protons to estimate their maximum energy. 

The essential point is that, if MFA to milligauss levels is occurring as part of DSA, 
the acceleration must be efficient and the system is strongly nonlinear. The accelerated 
particles and the pressure from the amplified field feedback on the shock structure (Fig. 14. Sp 
and this feedback makes the precursor less confining [i.e., {u{x)B[x)) <^ uqB2]- Therefore, 
a shock of a given physical size will not be able to accelerate protons to an energy as large 
as estimated ignoring NL effects. 

Despite the reduction in E^'^^ compared to test-particle predictions that our results 
imply, a remnant such as SNR RX J1713. 7-3946 might still produce CRs up to the knee. 
The NL example we have presented with B2 ~ 450 fxG produces protons up to ~ 100 TeV 
in ~ 100 yr in a confinement region of ~ 0.1 pc. If instead we had taken Lfeb = 1 pc, a size 
comparable to the western shell of SNR RX J1713. 7-3946, our NL model would produce 
~ IPeV protons in ~ 1000 yr. Protons of this energy are consistent with the ~ 30 TeV 
7-rays observed from SNR RX J1713. 7-3946 [68] and when the acceleration of heavy ions 
such as Fe^^^ is considered, the maximum particle energy extends to > 10^^ eV. 

As a final comment we emphasize a point also made by [21j . If MFA is occurring 
and the system is highly NL, it may not be possible to explain temporal variations in 
nonthermal X-ray emission simply as a radiation loss time. There cannot be variations 
in X-ray emission on short time scales unless the accelerator changes in some fashion on 
these time scales, otherwise the radiation would be steady, or varying on the shock dynamic 
timescale, regardless of how short the radiation loss time was. Since the injection and 
acceleration of protons and electrons is nonlinearly connected to the amplified magnetic 
field, changes in the electron particle distribution and changes in the field producing the 
synchrotron emission, will go together and it may be difficult to unambiguously determine 
the field strength from temporal variations. 
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4.3 Turbulence dissipation in shock precursor 

Here we present the results of our model regarding the effects of dissipation of 
turbulence upstream of the shock and the subsequent precursor plasma heatinjf. 

The magnetic turbulence generated by the instability is assumed to dissipate at a 
rate proportional to the turbulence generation rate, and the dissipated energy is pumped 
directly into the thermal particle pool (i.e., the model described by Equation (j3.73p is 
assumed). An iterative scheme is employed to ensure the conservation of mass, momentum, 
and energy fluxes, thus producing a self-consistent solution of a steady-state, plane shock, 
with particle injection and acceleration coupled to the bulk plasma flow modification and 
to the magnetic field amplification and damping. 

Our results show that even a small rate of turbulence dissipation can significantly 
increase the precursor temperature and that this, in turn, can increase the rate of injection of 
thermal particles. The nonlinear feedback of these changes on the shock structure, however, 
tend to cancel so that the spectrum of high energy particles is only modestly affected. 

4.3.1 Model 

We model the evolution of the turbulence, as it is being advected with the plasma 
and amplified, with the following equations: 

E±[U] = {l-aH)G±[U]+I±[U]. (4.7) 

Here, for readability, we abbreviated as E the evolution operator, as G the growth operator 
and as / the wave-wave interactions operator, acting on the spectrum of turbulence energy 
density U = {C/_(x, k), U+{x, k)}. These quantities are defined as follows: 

E±[U] = {u±Vg)-^U± + U±^(^u±Vg] , (4.8) 



dp 



dk 



(4.9) 



/±[C/] = ±^ ([/_-[/+). (4.10) 

The parameter an describes the turbulence dissipation rate, and for an = 0, the equa- 
tions (14. 7p becomes exactly the system of equations (14. ip and (14. 2p . In this system u = u{x) 



''Results presented here were a part of our article |123) 
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is the flow speed and Vq = Vg{x) is the parameter defining the turbulence growth rate and 
the wave speecj^. The parameter an enters the equations of turbulence evolution (j4.7p 
through the factor (1 — an), which represents the assumption that at all wavelengths only 
a fraction (1 — a^) of the instability growth rate goes into the magnetic turbulence, and 
the remaining fraction an is lost in the dissipation process. See Equation (j3.73p and the 
corresponding part of Section I3.2.4[ 

Equation (j3.76p is used to account for the precursor plasma heating by the dis- 
sipated turbulence. For L{x) = 0, equation (j3.76p reduces to the adiabatic heating law, 
-fth ~ p'^ and, for a non-zero L{x), it describes the heating of the thermal plasma in the 
shock precursor due to the dissipation of magnetic turbulence. 

The main effects of turbulence dissipation in our model are: (i) a decrease in the 
value of the amplified field Bcs{x), which determines the diffusion coefficient, D{x,p); (ii) 
an increase in the temperature of particles just upstream of the subshock, which infiuences 
the injection of particles into the acceleration process, and (iii) an increase in the thermal 
particle pressure Pth{x < 0), and a decrease in the turbulence pressure Pw{x), which enter 
the conservation equations described in Section 13.1.71 and 13.1.81 Since all of these processes 
are coupled, a change in dissipation influences the overall structure of the shock. 

4.3.2 Results 

Particle Injection in Unmodified Shocks (Subshock Modeling) 

In order to isolate the effects of plasma heating on particle injection, we first show 
results for unmodified shocks, i.e., u{x < 0) = uq and u{x > 0) = uo/'^tot, with fixed rtot- In 
these models particle acceleration, magnetic field amplification and turbulence damping are 
included consistently with each other, but we do not obtain fully self-consistent solutions 
conserving momentum and energy, since this requires the shock to be smoothed, while we 
intentionally fix u{x). 

In Fig. 14.101 we show results where the compression ratio is varied between rtot = 3 

and 3.6 as indicated. In all models, uq = 3000 km s^^, Tq = 10^ K, no = 0.3 cm~'^ and 

^As explained in [122) . in the quasi-linear case, AB <C -Bo, the wave speed and the speed determining 
turbulence growth rate are both equal to the Alfven speed, Vg{x) — va = Bq / \/ 4:Tt p{x) . In the case of 
strong turbulence, AB > Bo, we hypothesize that the resonant streaming instability can still be described 
by equations (|4.7p with Vc being a free parameter ranging from Bo/\/ 47rp(a;) to Bbs/ ^/4^Jvp{x). 




Figure 4.10: Dissipation effects in unmodified shocks 
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Bq = 3 fiG (the corresponding sonic and Alfven Mach numbers are Mso ~ Mao ~ 250). 
The FEB was set at xfeb = —3 • 10^ rgo (our spatial scale unit rgo = rripUQC / [eBo)) ^ and 
for each rtot we obtained results for different values of an between and 1. The values 
plotted in the top three panels of Fig. I4.1UI are the amplified magnetic field downstream, 
-Beff2) the Mach number right before the shock, M^i (this is not equal to M^q because of the 
plasma heating due to turbulence dissipation), and the fraction of thermal particles in the 
simulation that crossed the shock in the upstream direction at least once (i.e., got injected), 
/cr. The bottom panel shows the ratio of the calculated downstream effective magnetic field 
Bern to trend values -Btrend («//); what is meant by "trend" is the equation (j4.11j) explained 
in [123]: 



-Btrend(ai/) = f^O^'tot^) + (1 " "h) 



3/4^2 



(4.11) 



Looking at the curve for -Befr2 in the rtot = 3.0 and r^ot = 3.2 models, one sees 
an easy to explain behavior: as the magnetic turbulence dissipation rate, an-, increases, 
the value of the amplified magnetic field decreases, going down to B^r^^^ (the upstream 
field compressed at the shock) for uh = 1. Increasing an simply causes more energy to 
be removed from magnetic turbulence and put into thermal particles, thus decreasing the 
value of B^m ■ 

The plots for rtot ^ 3.4 present a qualitatively different behavior from those with 
rtot ^3.2. The downstream magnetic field -Bcfr2 does decrease with increasing a//, but not 
as rapidly as in the previous two cases, and there is a switching point at an ~ 0.95 in 
the curves for Mgi and /cr- The bottom panel of Fig. 14.101 shows a deviation of -Bcfr2 from 
the trend (j4.1ip by a large factor in the rtot = 3.4 case. This effect becomes even more 
dramatic for rtot = 3.5 and rtot = 3.6 where -Befr2) contrary to expectations, increases with 
an before an — > 1 . The fact that the final energy in turbulence can increase as more energy 
is transferred from the turbulence to heat indicates the nonlinear behavior of the system 
and shows how sensitive the acceleration is to precursor heating. 

It is worth mentioning that the observed increase of particle injection due to the 
precursor plasma heating is a consequence of the thermal leakage model of particle injection 
adopted here (see |123j for the explanation of this connection). In this model, a downstream 
particle, thermal or otherwise, with plasma frame speed v > U2, has a probability to return 
upstream which increases with v (see [9] for a discussion of the probability of returning 
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particles). An alternative model of injection (see, for example, [22]) is one where only 
particles with a gyroradius greater than the shock thickness can get injected. In the [22] 
model the fraction of injected particles may be insensitive to the precursor heating if the 
parameter controlling the injection in that model, ^, is fixed. While both models are 
highly simplified descriptions of the complex subshock (see, e.g., [901 [62]), they offer two 
scenarios for grasping a qualitatively correct behavior of a shock where particle injection 
and acceleration are coupled to turbulence generation and flow modification. Hopefully, 
a clearer view of particle injection by self-generated turbulence in a strongly magnetized 
subshock will become available when relevant full particle PIC or hybrid simulations are 
performed. 

With the general trends observed here in mind, we now show how nonlinear effects 
modify the effect dissipation has on injection and MFA. 

Fully Nonlinear Model 

In this section we demonstrate the results of the fully nonlinear models, in which 
the flow structure, compression ratio, magnetic turbulence, and particle distribution are all 
determined self-consistently, so that the fluxes of mass, momentum and energy are conserved 
across the shock. 

We use two sets of parameters, one with the far upstream gas temperature Tq = 
10^ K and the far upstream particle density no = 0.3 cm~'^, typical of the cold interstellar 
medium (ISM), and one with To = 10^ K and no = 0.003 cm~^ typical of the hot ISM. In 
both cases we assumed the shock speed uq = 5000 km s~^, and the initial magnetic field 
-Bo = 3 fiG (giving an equipartition of magnetic and thermal energy far upstream, nofc^To ~ 
Bq/^Stt)). The corresponding sonic and Alfven Mach numbers are Mg ~ Ma ~ 400 in both 
cases). The size of the shocks was limited by a FEB located at xfeb = — lO^rgo ~ —3 • 10~^ 
pc. For both cases, we ran seven simulations with different values of the dissipation rate 
an, namely oh G {0; 0.1; 0.25; 0.5; 0.75; 0.9; 1.0}. Also, for the hot ISM (To = lO'' K) case 
we ran a simulation neglecting the streaming instability effects, i.e., keeping the magnetic 
field constant throughout the shock and assuming that the precursor plasma is heated only 
by adiabatic compression (this model will be referred to as the 'no MFA case'). 

Tables HTT] and W?^ summarize some of the results of these models. The effect of the 
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Table 4.1: Summary of Non- linear Simulation in a Cold ISM 
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71 


21 
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See text and [123] for notation 



Table 4.2: Summary of Non-linear Simulation in a Hot ISM 
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See text and [123] for notation 
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turbulence dissipation into the thermal plasma is evident in the values of the pre-subshock 
temperature Ti, the downstream temperature T2, and the volume-averaged precursor tem- 
perature {T{x < 0)) (the averaging takes place between x = xfeb and x = 0). The value 
of Ti depends drastically on the level of the turbulence dissipation a^, increasing from 
an = to an = 0.5 by a factor of 100 in the cold ISM case, and by a factor of 11 in the 
hot ISM case. The values of the temperature as high as Ti are achieved upstream only 
near the subshock; the volume-averaged upstream temperature, {T{x < 0)), is significantly 
lower. The downstream temperature, T2, varies less with changing an, because it is largely 
determined by the compression at the subshock, which is controlled by many factors. It 
is worth mentioning the case without MFA reported in Table 14. 2[ Besides having a much 
larger compression factor than the shocks with MFA (r^ot = 13 as opposed to rtot ^8), it 
has a much smaller downstream temperature {T2 = 2.2-10^ K as opposed to T2 > 5.3-10^ K) 
These effects of dissipation on the precursor temperature may be observable. 

In Figure HTTT] we show results for /„, Mgi, -Befr2! and -Btrend which can be compared 
to the results for unmodified shocks shown in Figure I4.10[ For the modified shocks, the 
fraction of the thermal particles crossing the shock backwards for the first time, fcr, clearly 
increases by a large factor with an, which can be explained by the connection between 
Ti and the injection rate. One could expect that the amplified effective magnetic field 
BeS2 would behave similarly to the rtot = 3.5 case in Section 14.3.21 i.e. that -Beff2 would 
not decrease or even would increase for larger an- Instead, -Beff2 behaves approximately 
according to the trend ()4.1ip , as the values of -Btrcnd from Tables 14.11 and 14.21 show and the 
bottom panel of Fig. 14.111 illustrates. The important point is that, even though precursor 
heating causes the injection efficiency to increase substantially, the efficiency of particle 
acceleration (i.e., the fraction of energy in CRs) and magnetic turbulence generation is 
hardly changed. We base this assertion on the fact that i?efr2 remains close to -Btrend) which 
was derived under the assumption that changing an preserves the total energy generated 
by the instability, but re-distributes it between the turbulence and the thermal particles. 

Considering how much the injection rate fcr increases with an, and how much the 
upstream temperature of the thermal plasma, Ti, is affected by the heating, it is somewhat 
surprising that the trend of the amplified effective field B^m is unaffected. The mechanism 
by which the shock adjusts to the changing heating and injection in order to preserve the 
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MFA efficiency can be understood by looking at the trend of the total compression ratio 
rtot and the subshock compression ratio rgub in Tables 14.11 and 14. 2t they both decrease 
significantly for higher an- The decrease in rgub is easy to understand: with the turbulence 
dissipation operating in the precursor Mgi goes down, which lowers Tsui,. Additionally, 
decreasing P^i helps to reduce rgub, and with a boost of the particle injection rate, the 
particles returning for the first time increase in number and build up some extra pressure 
just upstream of the shock, which causes the flow to slow down in that region, thus reducing 
the ratio r^uh- 

Further understanding of the shock adjustment to the changing dissipation can 
be gained by studying Figures 14.121 - 14.151 in which we plot the spatial structure and the 
momentum-dependent quantities of the shocks in the cold ISM and the hot ISM cases for 
an G {0;0.5;1}. 

Figures [4 . 1 2 1 and 14 . 1 3 1 show an overlap in the curves for the flow speed u{x) in the 
an = and an = 0.5 models, differing only close to the subshock, where u{x) falls off more 
rapidly towards the subshock in the an = 0.5 case, resulting eventually in a lower r^ah- This 
means that for the high energy particles, which diffuse far upstream, the acceleration process 
will go on in about the same way with and without moderate turbulence dissipation (and 
the acceleration efficiency will be preserved with changing an)- For lower energy particles, 
however, there will be observable differences in the energy spectrum. The an = 1.0 case has 
a significantly smoother precursor, which is not unusual, given the lower maximal energy 
of the accelerated particles in this case (because of the magnetic field remaining low). The 
thermal gas temperatures T{x), plotted in the bottom panels of Figures [4. 121 and 14. 131 show 
that the temperature becomes high well in front of the subshock. 

The low energy parts of the particle distribution functions shown in Figures 14.141 
and 14.151 are significantly different for models with and without dissipation in both the 
cold ISM and the hot ISM cases. The apparent widening of the thermal peak refiects the 
increase in the downstream gas temperature T2. The differences extend from the thermal 
peak to mildly superthermal momenta 0.2 mpC, indicating an increased population of the 
'adolescent' particles with speeds up to w ~ 0.2c ~ 12uo when the turbulence dissipation 
operates. The high energy (p > 0.2 nipc) parts of the spectra for an = and an = 0.5 are 
similar (except for a lower Pmax due to a lower value of the amplified field in the an = 0.5 
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Figure 4.12: Nonlinear shocks with dissipation, cold ISM 
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Figure 4.13: Nonlinear shocks with dissipation, hot ISM 
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Figure 4.14: Particle distribution with dissipation, cold ISM 
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Figure 4.15: Particle distribution with dissipation, hot ISM 
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case), confirming our assertion about the preservation of the particle acceleration efficiency. 
The increased population of the low-energy particles just above the thermal peak should 
influence the shock's X-ray emission. 

The characteristic concave curvature of the particle spectra above the thermal 
peak is clearly seen in the top panels of Figures 14.141 and 14.151 These shocks are strongly 
nonlinear and, as the pressure spectra in the bottom panels show, most of the pressure is 
in the highest energy particles. For these examples, 60 to 80 percent of the downstream 
momentum flux is in CR particles. The number of particles producing this pressure is 
small, however, and as the plots in the middle panels show, the fraction of particles above 
the thermal peak is on the order of 10^'^, and the fraction of particles above 1 GeV is around 
10^^ in all cases. In addition to the pressure (and energy) in the distributions shown, a 
sizable fraction of shock ram kinetic energy flux escapes at the FEB. 

In Figure 14.161 the subshock region for the hot ISM case is shown enlarged for the 
models with (an = 0.5) and without dissipation (an = 0.0). For an = 0, the thermal pres- 
sure Pth remains low upstream (middle panel), and the subshock transition is dominated by 
the magnetic pressure Pw For an = 0.5 (bottom panel) the thermal pressure Pth just be- 
fore the shock becomes comparable with P^, but also the heating-boosted particle injection 
brings up the pressures of the 'adolescent' particles. For an = 0.5 the pressures produced 
by the first and second time returning particles (Pi and P2) are not small compared to Pth 
and Pw just upstream of the shock, which contributes to the reduction of rg^^ described 
above. However, the pressure of the 'mature' particles, P>5, doesn't change much, due to 
the non-linear response of the shock structure to the increased injection. 

To summarize, for both the unmodified (Fig. I4.10p and modified (Fig. 14. lip cases, 
Msi drops and fcr grows as an increases. The surprising result is that Pefr2 can increase 
in the unmodified shock as an goes up if rtot is large enough. This indicates that the 
boosted injection efficiency (i.e., larger fcr) outweighs the effects of field damping. This 
doesn't happen in the modified case (top panel of Fig. 14. lip because of the nonlinear effects 
from the increased injection. From Fig. 14.161 we see that the boosted injection results in 
a smoother subshock and this makes it harder for low energy adolescent particles to gain 
energy. Once particles reach a high enough momentum (p > 0.2mpc; see the top panel of 
Fig. I4.15P they diffuse far enough upstream where the boost in injection has a lesser effect. 
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Figure 4.16: Enlarged subshock region in the hot ISM case. 
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We must emphasize again that these results are very sensitive to the physics of 
particle injection at the subshock. It is difficult to predict how the nonlinear results would 
change if a different model of injection was used, but we can refer the reader to the analytic 
model [3] that uses the threshold injection model with a different diffusion coefficient. 

4.3.3 Discussion 

Our two most important results are, ffist, that even a small rate (~ 10%) of 
turbulence dissipation can drastically increase the precursor temperature, and second, that 
the precursor heating boosts particle injection into DSA by a large factor. The increase in 
particle injection modifies the low-energy part of the particle spectrum but, due to nonlinear 
feedback effects, does not significantly change the overall efficiency or the high energy part 
of the spectrum. Both the precursor heating and modified spectral shape that occur with 
dissipation may have observable consequences. 

The parameterization we use here is a simple one and a more advanced description 
of the turbulence damping may change our results. In our model the energy drained from 
the magnetic turbulence, at all wavelengths, is directly 'pumped' into the thermal particles. 
Superthermal particles only gain extra energy due to heating because the thermal particles 
were more likely to return upstream and get accelerated. In a more advanced model of 
dissipation, where energy cascades from large-scale turbulence harmonics to the short-scale 
ones, the low energy CRs might gain energy directly from the dissipation. It is conceivable 
that cascading effects might increase the overall acceleration efficiency, the magnetic field 
amplification, and the maximum particle energy a shock can produce. 

It is also possible that non-resonant turbulence instabilities play an important role 
in magnetic field amplification (e.g., [99|). This opens another possibility for the turbulence 
dissipation to produce an increase in the magnetic field amplification. For instance, [30] 
proposed a mechanism for generating long-wavelength perturbations of magnetic fields by 
low energy particles. If such a mechanism is responsible for generation of a significant 
fraction of the turbulence that confines the highest energy particles, then the increased 
particle injection due to the precursor heating may raise the maximum particle energy and, 
possibly, the value of the amplified magnetic field. 
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4.4 Bell's nonresonant instability and cascading in nonlinear 
model 

These results are currently in preparation for publication by Vladimirov, Bykov 
and Ellison. We will present the results of the model of nonlinear shock acceleration with 
amplification of strong stochastic magnetic fields by Bell's nonresonant streaming instabil- 
ity. We compare the assumption that the spectral energy transfer in the generated MHD 
turbulence is suppressed to the assumption that the Kolmogorov cascade determines the 
transfer. 

The results confirm that the nonresonant instability alone may produce a steady 
state shock structure with a very strong effective magnetic field. In addition, we find that, 
in the absence of cascading, the spectrum of the MHD turbulence is not a power law, as 
usually assumed, but has a prominent multiple-peak structure. The sharp peaks indicate the 
presence of eddies of different distinct scales. Also, the precursor of the shock is no longer 
smooth, but has several layers (i.e., it is stratified), where lower and lower energy cosmic rays 
are overtaken by the eddies and quickly accelerated. However, if the Kolmogorov cascade 
is assumed, the amplification of magnetic field is not as efficient, but the stratification is 
eliminated. 

We argue that the physically realistic solution is in between the two extreme cases 
that we presented here, and discuss the consequences of both scenarios for the process of 
particle acceleration by shocks and for the observable features of emitted radiation. 

4.4.1 Model 

We describe the evolution of turbulence by equation (j3.86p with boundary condi- 
tion (j3.87p . and set Ug = j3g = 0, and '^g = 5g = Sg = 1. For the turbulence amplification 
model, we choose the Bohm nonresonant instability, i.e., with T^^^ given by Equation (|3.69p . 

The flux of energy along the spectrum, H(x,A;), reflects the cascade of turbulent 
structures. Cascade of MHD turbulence may be anisotropic [64J, harmonics with wavenum- 
bers transverse to the uniform magnetic field experiencing a Kolmogorov-like cascade, while 
the cascade in wavenumbers parallel to the field is suppressed. The waves generated in the 
nonresonant instability are transverse, so the diffusion coefficient for particle transport par- 
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allel to the flow depends on the wavenumbers parahel to the magnetic field. It is uncertain 
whether the regime in which the instabihty operates wih lead to a Kolmogorov cascade, or 
to a suppression of the parallel cascade. We therefore consider two extreme cases: Model A, 
in which the cascading is fully suppressed, i.e., 11^ = 0, and Model B, in which the cascading 
is efficient and has the Kovazhny form (e.g., |120j ) given by equation (j3.79p . i.e., = Uk- 
The dissipation term, L, is assumed to be zero for Model A, and to have the form ()3.74p 
for Model B, i.e., Lb = Ly- For Model B, we also assume that the seed wave spectrum 
represents linear waves that are not subject to cascading or dissipation, and that the tran- 
sition to the turbulent regime takes place at a point xq where the amplified wave spectrum 
reaches the value kW{xQ, k) = Bq/At: at some k. At this point, 11 and L are set from zero to 
the values (|3.79p and (|3.74p . The wavenumber at which the dissipation begins to dominate, 
kd, is identified with the inverse of a thermal proton gyroradius: k^ = cBq / (cy^ rUpkBT) , 
where rrip is the proton mass, A;^ is the Boltzmann constant and T = T{x) is the local gas 
temperature determined from the gas heating induced by L, as described in Section [3. 2. 4[ 

Particle transport is described by the hybrid model of diffusion in strong turbu- 
lence, laid out in Section [331 

To calculate the diffusive current jd{x), we propagate the particles using the diffu- 
sion properties described above, and then compute the moment of the particle distribution 
function jd{x) = e f Vxf{x,p)d'^p by summing over all particles crossing certain positions. 

In order to determine the minimal particle gyroradius, rgi, that limits the long- 
wavelength generation by the instability as defined by Equation (j3.69p . we define the lowest 
CR momentum at the current position, pi, as the momentum below which the CRs con- 
tribute 1% of the total CR pressure. Then rgi is defined as rgi = cpi/{eBis) with Big 
calculated for the momentum pi . 

We use the iterative procedures (|3.58p and (j3.59p to achieve a self-consistent shock 
structure, in which particle distribution, turbulence spectrum and fiow structure are all 
consistent with each other, and the fundamental conservation laws are fulfilled. 
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4.4.2 Results 

We ran the Monte Carlo simulations of a shock with a speed uq = 10^ km s^^ 
propagating along a uniform magnetic field Bq = 3 fiG in a plasma with a proton density 
no = 0.3 cm^'^ and a temperature Tq = 10^ K. We assumed that the seed magnetic fluctua- 
tions have an effective value Ai?soed = Bq, and that the acceleration process is size- limited 
with a free escape boundary located at x = — lO"^ rgo, where rgo = muQc/eBo 3.5-10^° cm. 
Two simulations were performed, which as described in the previous section, we will call 
Model A and Model B. 
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Figure 4.17: Shocks with turbulence generation by Bell's nonresonant instability 
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The result for model A was the steady-state structure of a shock modified by 
efficient particle acceleration and magnetic field amplification, with a self-consistent com- 
pression ratio rtot = tio/'U2 ~ 15, a downstream magnetic field BQg{x > 0) 1000 fiG, and 
particle acceleration up to a maximum momentum Pmax ~ 10^ rUpC. Model B predicted 
a lower compression ratio, rtot ~ H, lower magnetic field Befi{x > 0) 120 fiG, and a 
maximum momentum Pmax ~ 2 • 10^ rUpC. 
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Figure 4.18: Spectral properties of shocks shown in Figure H.17I 

The self-consistent structure of the shocks (the flow speeds, the effective mag- 
netic field, the diffusive CR current and the thermal plasma temperature) are shown in 
Figure 14. 171 Besides the above mentioned differences in the compression ratio and the am- 
plified magnetic fields, one may notice the difference in the jd{x) plot. While the jdix) curve 
is smooth for Model B, it has an uneven structure for Model A, which reveals the stratifica- 
tion that becomes more apparent in the spectra of the self-generated magnetic turbulence 
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(see below). Another prominent difference is the significantly increased temperature T{x) in 
the precursor of the Model B shock, which comes about due to the dissipation of cascading 
turbulence at large k. 

In Figure [1.18t we show the particle distribution function f{p), the dependence of 
proton mean free path on momentum X{p), and the acceleration time to a certain momen- 
tum, t{p). The plots of f{p) show that shocks with either model of spectral energy transfer 
remain efficient particle accelerators: the concave shape indicates the nonlinear modifica- 
tion of the shock structure, also apparent in the plots of u{x). The thicker thermal peak 
and the higher low energy parts of the spectrum in the model with cascading are due to the 
increased turbulence dissipation, similarly to what was observed in |123j . The mean free 
path A(p) for model B is a smooth function of p, with X (x p for intermediate and X oc p'^ 
for the highest energy particles, but for model A it has plateaus that correspond to the 
trapping of particles by turbulent vortices of different scales. A similar uneven structure is 
seen in the acceleration time t{p): it has regions of rapid and slow acceleration. 
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Figure 4.19: Turbulence spectrum at different spatial locations. 
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The most intriguing result of this simulation is shown in Figure H.19[ Plotted there 
are the turbulence spectra, W{x, k), multiplied by k (so that a horizontal line represents a 
Bohm spectrum W oc k~^). As in the other figures, the dotted lines represent Model A, and 
the dashed lines - Model B, but here we also show the spectra at different locations, which 
adds lines of the same styles, but different thicknesses. The thickest lines correspond to the 
downstream region x > 0, the medium thickness lines - a point upstream of the subshock, 
X = —2 - 10^ Tgo, and the thinnest lines ~ to the unshocked interstellar medium (i.e., very far 
upstream). In both models, the spectra of stochastic magnetic fields are described by (|3.60p 
and represented by the thin horizontal lines. Closer to the shock, where a small current of 
streaming accelerated protons is present, fluctuations around k = 10^^ r^^ are amplified by 
the nonresonant instability. In model A, the energy spectrum of these fluctuations peaks 
around the value kc/2 corresponding to the maximum of Fnr from (I3.69p . but in model B, 
cascading spreads this energy over an extended inertial range of k (see the medium thickness 
lines). Closer to the shock, where lower energy particles appear, the generation of waves at 
k = 10~^ r~Q shuts down, according to the limits of applicability in (j3.69p . but the increased 
number of the streaming particles and accordingly raised diffusive current now corresponds 
to a greater kc, and shorter wavelength structures get amplified, around k = 10^^ r^Q. In 
model A, this results in a second peak of the turbulence spectrum at that wavenumber, 
but in model B, cascading smoothes out the spectrum. By the time the plasma reaches the 
downstream region (thick lines), three distinct peaks get generated with this mechanism in 
Model A, while Model obtains an amplified turbulence spectrum close to VF oc k^^. 

The peaks occur because of the coupling of particle transport with magnetic turbu- 
lence amplification. The first (smallest k) peak forms far upstream, where only the highest 
energy particles are present, and their current jd is low. These particles diffuse in the A oc 
regime, scattered by the short-scale magnetic field fluctuations that they themselves gen- 
erate. As the plasma moves toward the subshock, advecting the turbulence with it, lower 
energy particles appear. At some x, particles with energies low enough to resonate with 
the turbulence generated farther upstream (in the lowest k peak) become dominant. This 
strong resonant scattering leads to a high gradient oi jd (seen at x ~ —10^ r^Q in the third 
panel of Fig. I4.17p . and the wavenumber kc/2, at which the amplification rate Fnr has a 
maximum, increases rapidly. The increased value of kc/2 leads to the emergence of the 
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second peak between 10~^ and 10~^ r^^, as seen in Fig. I4.19[ Similarly, the third peak is 
generated at distances closer to the subshock than ~ —2 x 10^ rgo and this is seen in the 
thick dotted line in Fig. 14.191 at k ~ 10r~o^ 

The number of peaks depends on the dynamic range, i.e., on Lfeb- A smaller 
LpEB can result in two peaks, while a larger Lfeb, and therefore a larger Pmax, can yield 
four or more peaks in the downstream region. 

The formation of the spectrum with discrete peaks occurs simultaneously with the 
stratification of the shock precursor into layers (see the plots of j^i), in which vortices of 
different scales are formed. The peaks are a direct result of Bell's nonresonant instability, 
but they will not show up unless Ass and Fnr are calculated consistently, and the simulation 
has a large enough dynamic range in both k and p. 

4.4.3 Discussion 

Our results show that, similarly to the model with a Bohm diffusion coefficient and 
a resonant streaming instability |122l 1123] . the predictions of efficient particle acceleration, 
shock modification and magnetic field amplification by a large factor remain in force. How- 
ever, compared to the previous results, the precursor structure is strikingly different: instead 
of a smooth, gradual variation of all quantities in the precursor, we observe a stratification 
into layers, in which vortices of distinct sizes are subsequently generated. The resulting 
turbulence spectrum has 3 sharp peaks, including one at very short wavelengths. This 
stratification process is eliminated if the rapid Kolmogorov cascade of turbulence structures 
is assumed. In the latter case, the amplified turbulence spectrum downstream becomes a 
power law W tx k~^, and the variation of all quantities in the shock precursor reverts to 
being smooth. The amplified effective magnetic field, the shock compression ratio and the 
maximum energy of accelerated particles are smaller in the model with Kolmogorov cas- 
cade. Cascading of MHD turbulence with respect to the wavenumbers parallel to the mean 
magnetic field must be suppressed [64], and the two models (without cascading and with 
the Kolmogorov cascade) should be perceived as the extremes, between which the more 
physically realistic answer lies. 

If the situation without cascades and with precursor stratification is the better ap- 
proximation of physical reality, what consequences for astrophysical observations and theory 
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of cosmic accelerators might it have? The calculation performed here derived a steady-state 
structure of a size-limited particle accelerator, but the information about mean particle ac- 
celeration time, Tacc(p) (bottom panel of Figure [d8j) allows a peek into the time-dependent 
process. The time of acceleration to a certain momentum is, on the average, proportional 
to the momentum, Tacc oc p, for p > nipC is, but there are periods of slow acceleration, 
when dluTacc/dlnp > 1, and fast acceleration, when dlnracc/rflnp < 1. Does it mean that 
in a time-dependent calculation, one would observe quiet periods, when the highest energy 
particles escape ahead of the shock into the interstellar medium and generate large-scale 
turbulent vortices, intermittent with bursts of particle acceleration, when the lower energy 
particles are trapped by these vortices and vigorously accelerated? The large amount of en- 
ergy observed in the shortest-scale (largest k) peak may influence the synchrotron radiation 
of electrons. Its location was around 0.1 Vgo ~ 3.5 • 10^ cm, and it contained roughly 1/3 of 
the magnetic field energy corresponding to the 1000 /iG magnetic field. Will this rapidly 
varying field affect the radio or the X-ray (e.g., [31j ) part of the SNR shock synchrotron 
spectrum? 

On the other hand, if the Kolmogorov cascade is the better representation of the 
spectral energy transfer in the problem of diffusive shock acceleration, it means that the 
amplified magnetic fields may not be as large as the quasi-linear theory suggests. Also, the 
heating of the shock precursor by the dissipation of turbulence must have significant effects. 
Indeed, in the model with cascades, the upstream plasma temperature T(x) is increased to 
values above 10^ K (bottom panel of Figure [4.17p . and the accelerated particle spectrum 
is elevated up to p ~ nipC with respect to model A. These features of the solution indicate 
that the X-ray emission of shocks must carry the fingerprints of the turbulent cascade. 
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4.5 Fits for the nonlinear shock structure 

The predictions of the model: how efficient magnetic field amplification and par- 
ticle acceleration are, how much the nonlinearly modified shocks compress and heat the 
medium that they propagate in, are important for many applications where strong shocks 
exist. Although nonlinear shock structure is likely to emerge in many problems, there exists 
no simple description for physicists working in other areas for to incorporating the nonlinear 
effects into calculations. 

In an effort to fix this situation, I derived simple scaling laws to replace the Hugo- 
niot adiabat, when nonlinearly modified shocks are considered (see also one such scaling in 
[25j). I did it by performing the derivation of the self-consistent shock structure for a set 
of input parameters spanning a certain range. After that, using the least squares method, 
I derived the best fit coefficients for power-law scalings fitting the obtained data. 

4.5.1 Model 

The nonlinear shock model used here is identical to that described in Section [4.31 

4.5.2 Results 

I ran a total of 81 Monte Carlo simulations with different parameters and obtained 
a self-consistent solution in each case. I chose the parameter range that represents the 
conditions in galaxy gluster shocks: an = . . . 1, Tq = 2 x 10'^ K, Bq = 0.1 ... 1.0 fiG, 
no = 10~^ . . . 10~^ cm~'^, and uq = 1000 . . . 3000 km s^^. For densities at the lower end of 
the range, magnetic field was only varied between 0.1 /xG and 0.5 — 0.6 ^G, because Alfven 
Mach numbers are too low for high Bq = 1 /iG for low density. The upstream free escape 
boundary was chosen as xfeb = — 10''rgo. 

The raw data I collected are shown in the tables IT3l 14.41 and 14.51 In these tables, 
the input parameters of the models are listed to the left of the vertical divider: an is the 
dissipation rate parameter, uq is the shock speed, Tq is the far upstream gas temperature, no 
is the far upstream plasma density, and Bq is the far upstream magnetic field. The rest of the 
columns are the self-consistent results of the simulation: rtot and rgub are the self-consistent 
total and subshock compression ratios, Ti is the temperature right before the subshock. 
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Bern is the downstream amplified effective magnetic field and T2 is the downstream gas 
temperature. 

I fitted T2 and rtot from these data with power law fits in the form: 

r2 = ^<^i?o"2<3^ (4.12) 
rtot = Cu^'B^o^n^Q^ (4.13) 

I chose to fit the an = 0, an = 0.5 and an = 1 cases separately. The results are shown 
below. The temperature scales as 



T2 

10«K 
106 K 



X 1.40±0.17 / o X 0.47±0.10 . -0.22±0.08 

i.n ( _^y-™ /^y-" - -0.3.0.03 

X 1.40±0.06 / ^ X 0.44±0.03 „^ s -0.21±0.03 



The compression ratio scales as 



x0.34±0.10 / 5 X -0.25±0.06 „^ .0.12±0.05 



„^ X 0.37±0.06 / N -0.25±0.03 „^ .0.13±0.03 



= 12.0 ^^^;3^j ^—J (_j , (4.18) 

/ «o \0.35±0.04 . ^ X -0.25±0.02 0.12±0.02 ^ ^ 

The deviations of the parameters shown above come from the standard least squares method 
and are 2a (95% confidence). To estimate the fit quality, I also calculated the mean square 
relative error and the maximum error of the fits in each case. For the temperature fits the 
mean square deviations of the fits (j4.14p . (j4.15p and (j4.16p from the data were 18%, 7% and 
6%, respectively, and the maximum errors were 41%, 25% and 17%, respectively. For the 
compression ratio fits the mean square deviations of the fits (j4.17p . (j4.18p and (j4.19p from 
the data were 12%, 7% and 5%, respectively, and the maximum errors were 42%, 26% and 
%, respectively. 



4.5.3 Discussion 

I calculated the self-consistent structure of nonlinear shocks that power particle 
acceleration and magnetic field amplification. The parameter range I spanned makes these 
calculations applicable to cosmological shocks [26], except for free escape boundary location, 
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^^FEB- Galaxy cluster formation shocks may have much larger spatial scale than defined 
by xpEB, but running the Monte Carlo simulation with a much greater xfeb is too time 
consuming. However, I argue and cans show with simulation results that as soon as xfeb is 
large enough to ensure that the fraction of ultra-relativistic particles in the shock precursor 
is significant, increasing xfeb further does not affect the self-consistent compression ratio 
and the downstream temperature too much (see also [123j ). 

By fitting the results of a number of simulations, I derived simple scaling laws for 
the downstream temperature and the shock compression ratio, expressed by the equations 

([4141) - dun]). 

These predictions are, of course, very different from the hydrodynamic shock so- 
lution. For instance, the Hugoniot adiabat ()3.39p and (|3.40p provides the following scaling 
for Ms > 1: 

(-^) = 1.2- 10^ f X , (4.20) 

V106Ky VlO^kms^V 

rtot = 4. (4.21) 

In the fits that I found, the temperature is orders of magnitude lower, and the compression 
ratio several times greater, than in the hydrodynamic shock solution. 

This method can be extended to different parameter ranges, and our model can 
make similar predictions of macroscopic parameter scalings for shocks in other systems. For 
example, the emission spectra of radiative shocks (from the infrared to the X-ray ranges) 
may be influenced by particle acceleration and magnetic field amplification. 

I am grateful to A. M. Bykov for the idea of this direction of research. 
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Table 4.3: Self-consistent shock parameters for an = 0.0 



an 


uo, km s 




To, K 


no 


, cm 


Bo, G 


r-tot 


^sub 




Ti, K 


Bc«2, G 




T2, G 


0.0 


1 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-07 


7 


19 


2 


87 


3 


63EH 


-04 


2 


25E-06 


2 


70E+06 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-07 


8 


50 


2 


91 


4 


06EH 


-04 


3 


84E-06 


7 


94E+06 


0.0 


3 


OE+03 


2 


OE+0-1 


1 


OE-04 


1 


0E-07 


8 


96 


2 


94 


4 


25EH 


-04 


5 


13E-06 


1 


64E+07 


0.0 


1 


OE-l-03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


6 


61 


2 


86 


3 


50E+04 


2 


76E-06 


3 


16E+06 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


7 


13 


2 


89 


3 


71E+04 


4 


70E-06 


1 


lOE+07 


0.0 


3 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


8 


12 


2 


92 


4 


09E+04 


6 


82E-06 


1 


94E+07 


0.0 


1 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-06 


3 


92 


2 


81 


2 


54E+04 


3 


30E-06 


8 


29E+06 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


5 


69 


2 


89 


3 


28E+04 


5 


74E-06 


1 


70E+07 


0.0 


3 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


6 


39 


2 


91 


3 


69E+04 


8 


20E-06 


3 


07E+07 


0.0 


1 


OE+03 


2 


OE+04 


3 


0E-05 


1 


OE-07 


6 


95 


2 


85 


3 


60E+04 


1 


31E-06 


2 


83E+06 


0.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


8 


05 


2 


92 


4 


OOE+04 


2 


57E-06 


8 


90E+06 


0.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


9 


01 


2 


91 


4 


37EH 


-04 


3 


42E-06 


1 


57E+07 


0.0 


1 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


5 


68 


2 


83 


3 


22EH 


-04 


1 


65E-06 


4 


13E+06 


0.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


6 


81 


2 


90 


3 


66EH 


-04 


2 


98E-06 


1 


20E+07 


0.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


7 


20 


2 


92 


3 


90EH 


-04 


4 


14E-06 


2 


43E+07 


0.0 


1 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


3 


62 


2 


80 


2 


42EH 


-04 


1 


85E-06 


9 


53E+06 


0.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


5 


60 


2 


88 


3 


28EH 


-04 


3 


19E-06 


1 


73E+07 


0.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


6 


13 


2 


92 


3 


67EH 


-04 


4 


52E-06 


3 


32E+07 


0.0 


1 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


6 


57 


2 


83 


3 


56EH 


-04 


8 


83E-07 


3 


llE+06 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


7 


52 


2 


90 


3 


86EH 


-04 


1 


60E-06 


9 


99E+06 


0.0 


3 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


8 


00 


2 


91 


4 


09EH 


-04 


2 


24E-06 


1 


99E+07 


0.0 


1 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


4 


02 


2 


82 


2 


59EH 


-04 


1 


02E-06 


7 


89E+06 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


5 


96 


2 


88 


3 


41E+04 


1 


83E-06 


1 


55E+07 


0.0 


3 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


6 


41 


2 


91 


3 


77E+04 


2 


54E-06 


3 


03E+07 


0.0 


1 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


2 


63 


2 


60 


2 


OlE+04 


1 


20E-06 


1 


47E+07 


0.0 


2 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


4 


51 


2 


88 


2 


95E+04 


1 


95E-06 


2 


62E+07 


0.0 


3 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


5 


49 


2 


92 


3 


51E+04 


2 


75E-06 


4 


09E+07 
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Table 4.4: Self-consistent shock parameters for an = 0.5 



OiH 


ito, km s 




To, K 


no 


, cm 


Bo, G 


r-tot 


^sub 




Ti, K 


Bc«2, G 




T2, G 


0.5 


1 


OE+03 


2 


OE+04 


1 


0E-04 


1 


OE-07 


6.68 


2 


36 


1 


28EH 


-06 


1 


54E-06 


3 


30E+06 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-07 


7.95 


2 


55 


1 


95EH 


-06 


2 


36E-06 


8 


68E+06 


0.5 


3 


OE+03 


2 


OE+0-1 


1 


0E-04 


1 


OE-07 


8.74 


2 


62 


2 


87EH 


-06 


3 


llE-06 


1 


62E+07 


0.5 


1 


OE-l-03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


5.39 


2 


24 


2 


18E+06 


1 


68E-06 


4 


91E+06 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


6.73 


2 


34 


4 


71E+06 


3 


07E-06 


1 


24E+07 


0.5 


3 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


7.44 


2 


42 


7 


49E+06 


4 


43E-06 


2 


27E+07 


0.5 


1 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


3.65 


2 


40 


2 


40E+06 


2 


99E-06 


8 


58E+06 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


4.96 


2 


26 


9 


29E+06 


4 


19E-06 


2 


18E+07 


0.5 


3 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


5.72 


2 


26 


1 


67E+07 


5 


66E-06 


3 


82E+07 


0.5 


1 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


5.98 


2 


29 


1 


68E+06 


8 


36E-07 


3 


97E+06 


0.5 


2 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


7.45 


2 


40 


3 


53E+06 


1 


59E-06 


1 


02E+07 


0.5 


3 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


8.32 


2 


48 


4 


85EH 


-06 


2 


lOE-06 


1 


77E+07 


0.5 


1 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


4.66 


2 


21 


2 


74EH 


-06 


1 


16E-06 


6 


06E+06 


0.5 


2 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


5.94 


2 


25 


7 


20EH 


-06 


1 


92E-06 


1 


59E+07 


0.5 


3 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


6.63 


2 


29 


1 


21EH 


-07 


2 


77E-06 


2 


89E+07 


0.5 


1 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


3.44 


2 


46 


2 


04EH 


-06 


1 


72E-06 


9 


52E+06 


0.5 


2 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


4.82 


2 


25 


9 


95EH 


-06 


2 


43E-06 


2 


28E+07 


0.5 


3 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


5.58 


2 


30 


1 


68EH 


-07 


3 


14E-06 


3 


98E+07 


0.5 


1 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


5.36 


2 


20 


2 


30EH 


-06 


5 


46E-07 


5 


05E+06 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


6.70 


2 


29 


5 


58EH 


-06 


1 


04E-06 


1 


29E+07 


0.5 


3 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


7.52 


2 


36 


7 


68EH 


-06 


1 


42E-06 


2 


16E+07 


0.5 


1 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


3.78 


2 


38 


2 


47EH 


-06 


9 


24E-07 


8 


21E+06 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


4.96 


2 


27 


9 


51E+06 


1 


30E-06 


2 


19E+07 


0.5 


3 


OE+03 


2 


OE+04 


1 


OE-05 


3 


OE-07 


5.77 


2 


26 


1 


64E+07 


1 


77E-06 


3 


81E+07 


0.5 


1 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


2.62 


2 


58 


1 


41E+05 


1 


20E-06 


1 


46E+07 


0.5 


2 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


4.17 


2 


37 


9 


86E+06 


1 


68E-06 


2 


84E+07 


0.5 


3 


OE+03 


2 


OE+04 


1 


OE-05 


5 


OE-07 


4.93 


2 


30 


2 


04E+07 


2 


09E-06 


4 


86E+07 
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Table 4.5: Self-consistent shock parameters for an = 1.0 



an 


uo, km s 




To, K 


no 


, cm 


Bo, G 


r-tot 


^sub 




Ti, K 


Bc«2, G 




T2, G 


1.0 


1 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-07 


5.86 


2 


12 


2 


63EH 


-06 


3 


82E-07 


4 


72E+06 


1.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-07 


7.52 


2 


29 


4 


66EH 


-06 


4 


57E-07 


1 


04E+07 


1.0 


3 


OE+03 


2 


OE+0-1 


1 


OE-04 


1 


OE-07 


8.35 


2 


44 


5 


96EH 


-06 


4 


93E-07 


1 


84E+07 


1.0 


1 


OE-l-03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


5.08 


2 


07 


3 


46E+06 


1 


03E-06 


5 


88E+06 


1.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


5.91 


2 


02 


1 


lOE+07 


1 


16E-06 


1 


82E+07 


1.0 


3 


OE+03 


2 


OE+04 


1 


OE-04 


3 


OE-07 


6.96 


2 


23 


1 


36E+07 


1 


30E-06 


2 


80E+07 


1.0 


1 


OE+03 


2 


OE+04 


1 


OE-04 


1 


0E-06 


3.43 


2 


25 


4 


06E+06 


2 


74E-06 


1 


OOE+07 


1.0 


2 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


4.49 


2 


07 


1 


66E+07 


3 


18E-06 


2 


89E+07 


1.0 


3 


OE+03 


2 


OE+04 


1 


OE-04 


1 


OE-06 


5.17 


2 


09 


3 


OlE+07 


3 


50E-06 


5 


12E+07 


1.0 


1 


OE+03 


2 


OE+04 


3 


0E-05 


1 


OE-07 


5.48 


2 


07 


3 


07E+06 


3 


63E-07 


5 


16E+06 


1.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


6.87 


2 


25 


5 


70E+06 


4 


26E-07 


1 


24E+07 


1.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


1 


OE-07 


7.69 


2 


28 


9 


80EH 


-06 


4 


67E-07 


2 


22E+07 


1.0 


1 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


4.25 


2 


08 


4 


74EH 


-06 


9 


19E-07 


7 


98E+06 


1.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


5.32 


2 


08 


1 


30EH 


-07 


1 


07E-06 


2 


19E+07 


1.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


3 


OE-07 


5.96 


2 


10 


2 


18EH 


-07 


1 


16E-06 


3 


89E+07 


1.0 


1 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


3.32 


2 


38 


3 


47EH 


-06 


1 


63E-06 


1 


05E+07 


1.0 


2 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


4.42 


2 


12 


1 


67EH 


-07 


1 


90E-06 


2 


93E+07 


1.0 


3 


OE+03 


2 


OE+04 


3 


OE-05 


6 


OE-07 


5.01 


2 


11 


3 


21EH 


-07 


2 


06E-06 


5 


47E+07 


1.0 


1 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


4.95 


2 


04 


3 


85EH 


-06 


3 


40E-07 


6 


33E+06 


1.0 


2 


OE+03 


2 


OE+04 


1 


OE-05 


1 


OE-07 


6.04 


2 


12 


1 


OlEH 


-07 


3 


92E-07 


1 


75E+07 


1.0 


3 


OE+03 


2 


OE+04 


1 


OE-05 


1 
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4.6 Spectrum and angular distribution of escaping particles 

Using our model, we calculated the spectra of particles escaping from the shock 
at the free escape boundary. We provide simple fits to the energy and angular distribution 
of the escaping particles. 

4.6.1 Model 

For Bohm diffusion, the momentum distribution f(p) of the escaping particles can 
be described as 



P/Pn 



/ f dx/x 

exp —s j — '— 



m oc p-' \ \ , ' . (4.22) 

exp (Pmax/P) - 1 

Here s is the power-law index corresponding to the compression ratio r, s = 3r/(r — 1) 
and Pmax is defined below. The numerator describes the exponential turn-over of the high 
energy particles, and the denominator describes the low-energy part of the escaping particle 
distribution. Equation (I4.22|) is similar to equations (7) and (8) of |129j . but assumes a 
diffusion coefficient D{p) oc p, as opposed to D{p) oc assumed in [129] . 

To find the normalization of the escaping particle distribution f{p), one needs to 
use the quantity gesc self-consistently defined by the simulation as ()3.53p : 

oo 

47r / p'^dp / dfj,f{p)g{fi)cp = -QescPoul, (4.23) 



-1 

where (/esc is the fraction of energy flux carried away by escaping particles, and g{fi) is their 
angular distribution. The latter function is defined so that J^^ g{fi)dfi = 1, and fi = Px/p- 
The quantity pmax (the maximum particle momentum) can be estimated from the 
test-particle theory of particle acceleration as Pmax ~ 3noei?o|2;FEB|/c^ (i-e., the momentum 
at which the Bohm diffusion length equals |xfeb|)- The function g(fj,) is the distribution of 
particles incident on a fully absorbing boundary in a flow moving at a speed u. It can be 
estimated using the Monte Carlo simulation, as shown below. 
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4.6.2 Results 

I ran a simulation of a nonlinearly modified shock with a speed uq = 5000 km s^^, 
compression ratio r ~ 10, no magnetic field amplification, and Bohm model of diffusion. 
Then I plotted and fitted the particle distribution (angular and momentum-space) deter- 
mined by the simulation at the free escape boundary located at xfeB) as as shown in 
Figures 14.201 and 14.211 The histograms shown in these figures are the results of the Monte 
Carlo simulation, and the smooth lines are the fits, equations for which are provided in the 
figures. 




Figure 4.20: Angular distribution of escaping particles 



p"'/'f(p) ~ exp(-(10/3)\inU0|°|p/p„|j1/({e'^''-1)x)dx!) / {exp(pyp)-l) 
p„ = 1 .2 X 3u„eBoXpEB/c^ 
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Figure 4.21: Momentum distribution of escaping particles 



It turns out that the angular distribution of the particles can successfully be fitted 
with the equation shown in Figure 14.201 



0.70|/i|2 -h0.65|/i|, if /X < 0, 
0, if /i > 0. 



(4.24) 
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This represents the case when all the particles are moving to the left, i.e., away from the 
shock (because g{fx > 0) = 0). The momentum distribution of escaping particles has the 
shape described by equation (j4.22p with 

SuoeBo 

Pmax = 1-2 5 \xfEb\- (4.25) 

The factor 1.2 is a minor correction to the above mentioned analytic estimate. 

I would like to thank T. Kamae and S.-H. (Herman) Lee for the discussions that 
have led to the results presented here. 

4.6.3 Discussion 

I used our model to quantitatively describe the spectrum and the anisotropic 
angular distribution of particles leaving the shock at the upstream free escape boundary. 
The calculations accounted for the nonlinear modification of the shock by efficient particle 
acceleration. 

While the particular results derived here have limited applicability, because they 
apply to just a single set of parameters of a plane shock, I provided them to indicate a 
possible direction of research applying the model presented in this dissertation. 

One application may be the description of the interaction between the CRs escap- 
ing from a shock and the interstellar medium. For instance, the streaming of these particles, 
carrying a large fraction of the shock's energy flux, may amplify magnetic field fluctuations 
in the interstellar medium. 

Another interesting astrophysical application of these results is the recently dis- 
covered diffuse gamma ray sources identified as molecular clouds illuminated by cosmic ray 
particles produced in a nearby supernova remnant shock wave (e.g., [88j). A major un- 
certainty for the interpretation of these observations is whether the particle accelerating 
shocks are traversing through the cloud, or located far away from it (see, e.g., [59]). In 
the first case, the angular distribution of accelerated protons is far from isotropic, while in 
the second case the CRs may have had time to isotropize in the interstellar medium before 
they reach the cloud. The gamma ray emission of these protons via the decay of pi-mesons 
produced in collisions with the cloud gas protons will be different in these two cases, because 
the relativistic process produces gamma rays with strong dependence of energy spectrum 
on the angle of emission. 
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Chapter 5 

Conclusions 

I have developed a model of nonlinear shock acceleration that self-consistently 
includes the amplification of stochastic magnetic fields in the shock precursor by the accel- 
erated particles produced in the first order Fermi process. The model is based on the Monte 
Carlo simulation of particle transport developed by Ellison and colleagues, and my contribu- 
tion to the model was the incorporation of the analytic description of magnetic turbulence 
amplification and evolution, and the implementation of particle transport consistent with 
the generated magnetic turbulence. 

In this dissertation, I provided the details of the model to a degree that, I believe, 
make it reproducible. I presented the tests of the various parts of the simulation, which 
confirm that the results of the computer code I built agree with the known analytic results. 
This dissertation also contains an outline of our three refereed publications, in which we 
presented the applications of our model. It also features some results that have not yet 
been published. 

The applicability of the Monte Carlo model to shock acceleration in space was 
tested well before my work in this project by Ellison, Baring and others [43, 8j, who used 
the most direct data available - spacecraft observations of the Earth's bow shock and of 
interplanetary shocks. The physical correctness of the magnetic field amplification that I 
implemented is yet to be tested. Nevertheless, the predictions of the model are able to 
explain the observations that inspired it (i.e., the large magnetic fields and increased shock 
compression ratios in SNRs). 
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The most important limitations of the model are the uncertainty of the extrapola- 
tions of linear models of magnetic turbulence evolution into the nonlinear regime, and the 
statistical description of particle transport in stochastic magnetic fields, which is subject to 
various conjectures. However, the strength of the presented model compared to the simpli- 
fied analytic treatments of particle acceleration and magnetic field amplification in shocks is 
the self-consistency. Our model allows one to determine the shock structure, the accelerated 
particle spectrum and the turbulence generation all consistently with each other, through an 
iterative procedure. Even compared to the advanced analytic nonlinear models of particle 
acceleration (e.g., [3]), our simulation stands out because the Monte Carlo technique can 
handle anisotropic particle distributions, which is essential for a more precise description of 
plasma physics; for instance, the injection of particles into the acceleration process is pre- 
dicted self-consistently in our model, yet it requires an additional parameter in others. Also, 
the inclusion of various factors that determine the plasma physics (e.g., turbulent cascades) 
is straightforward in our approach, and may be complicated in the analytic calculations. 

The applications of the model developed here are numerous and exist wherever 
strong non-relativistic collisionless shocks are present. This includes, to some degree, shocks 
in interplanetary space, supernova remnants, shocks in galaxy clusters, etc. Our results help 
answer questions regarding the sources of galactic cosmic rays up to the 'knee' of the CR 
spectrum, and they may be used in the modeling of supernova remnants, galaxy clusters 
and other objects. 

The dissertation contains the results of the model in Chapter HI Sections 14.11 
presents our first refereed publication [122j featuring the model. In this work we studied 
the self-consistent structure of particle accelerating shocks in the presence of the resonant 
CR streaming instability (see Section 13.2. ip . We confirmed that the efficient particle accel- 
eration and strong magnetic field amplifications can exist in collisionless shocks in a wide 
range of the possible rates of nonlinear development of the streaming instability. In Sec- 
tion [321 I present our article [55j that discusses the impact of magnetic field amplification 
on the maximum energy of the accelerated particles. We showed that the amplified mag- 
netic field does increase the maximum achievable particle energy, but by a smaller factor 
the increase of the field. Section [4.31 contains the results of our investigation of the effect of 
turbulence dissipation in the shock precursor. These results, presented in [123], show that 



173 



the conversion of turbulent energy into heat increases the pre-shock temperature, which af- 
fects particle injection into the acceleration process. In addition to presenting the published 
articles, I included some work in progress in this dissertation. In Section [4.41 1 demonstrate 
the results of the simulation of the nonlinear shock structure with the nonresonant Bell's 
instability (see Section I3.2.2p and the hybrid model of particle diffusion (see Section I3.3.5P . 
We find that, if turbulent cascade is suppressed, the self-consistent steady state shock struc- 
ture has a stratified precursor, and the spectrum of turbulence has an unusual multiple-peak 
structure. Section [4. 51 contains an outline of the calculations that can be done with the sim- 
ulation in order to obtain simple power-law fits to the results of nonlinear DSA. Such fits 
can be used in the models of supernova remnants, galaxy cluster shocks, or other objects 
where strong particle-accelerating shocks are present. Finally, in Section [4.61 1 describe how 
the escaping particle distribution can be fitted with simple functions and where these fits 
may be used. 

Considering the rapid growth of observational X-ray and gamma ray facilities that 
reveal the 'high energy Universe', such as Chandra, Fermi, H.E.S.S., etc., I believe that the 
development of this model is very timely and beneficial for the research in different areas 
of astrophysics. 
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Appendix A 

Numerical integrator for model with isotropization 

In our works [122j , [53] and [123J , we were including the generation by the stream- 
ing instability of waves traveling in both direction, and used the version of the wave ampli- 
fication equation that accounts for the interaction between these waves. In fact, this effect 
is only important for the weaker shocks that we did not consider, and in the more recent 
version we neglected the waves traveling downstream. In order to make a record of the 
previous work, I provide here the numerical integrator of the previously used model. 

The equations we will now consider are (j4.1|] and (14. 2p . They include the resonant 
streaming instability (generating and damping waves traveling in both direction) , the effect 
of wave amplitude increase with plasma compression, and the nonlinear interactions between 
the waves traveling in different directions. 

These equations, by introducing quantities 



C{x,kj) = J {U-{x,k) + U+{x,k))dk, (A-1) 
rj{x,kj) = [ {U^{x,k) -U+{x,k))dk, (A-2) 



are transformed into 

u^' -Vwv' + -v'^T] -jv^,^ = 0, (A-3) 

3 dP 2 

uv' -Vw^' + i:u'v-v'^C-v^— + —rj = 0. (A-4) 

2 dx Tr 

Here dP/dx is the gradient of CR pressure produced by particles resonant with waves in 
the bin Afcj (see resonance condition below). Then the latter are re- written as 

Aif + By + c = Q, (A-5) 

where 



u —Vu, . 

A= \ \ , B 

-Vw U 





(A-6) 
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If equation ()A-5P was linear (i.e., the matrices A, B and c did not depend on y), 
then solving it would be straightforward. I will skip the details and leave it to the reader 
to verify that the solution of equation (|A-5|) with the initial condition 

y(0) = yo, (A-7) 
assuming that A is reversible (otherwise, the system is not consistent) is 

yo{x) = exp [-A'^Bx] yo - 

However, c explicitly depends on y, and for fwf > 0, the magnetic field determining Vu, 
depends on the integral of U± with respect to k (see Section [4.ip . making the matrices A 
and B depend on y in a non-trivial way, and therefore a numerical solution is required. 

To integrate ()A-5p . let us start off by assuming that y{x) = yo for any x. Then 
let the integrator perform a ^level-f iterative procedure, the purpose of which is to deal 
with the fact that in (IA-5|) the matrices A and B depend on the unknown functions (,{x, k), 
rj{x,k) in all x-space and fc-space through Vw{x) depending on U±{x,k). Here is what 
the Hevel-V iterative procedure involves. Given a fc-bin, integrate the equations (jA-Sh for 
that bin from far upstream to downstream. The values of U-{x,k) and U-^-{x,k) used to 
form matrices A and B are the ones obtained from the previous iteration. After all /c-bins 
have been integrated, the iteration is over. Then the just obtained values of y{x, k) = 
{S,{x,k), ri{x,k))'^ are used to run the next iteration. Iterating on Hevel-V ends when the 
current iteration gives results that are close enough to the results of the previous iteration. 

Integrating from one grid plane (at xq) to the next one (at xi), the routine is not 
likely to encounter a very strong variation of determining the matrices A and B (because 
the latter depend on an integral of y with respect to k) , but the quantities r/ and ^ that enter 
c may vary by a large factor, and care must be take with using (jA-Sp . I employ another 
iterative procedure {Hevel-2'') to tend to the dependence of the vector c, on (^{x,k), r]{x,k) 
in (|A-5p . This iterative procedure is described below. First, divide the step from xq to xi 
into Nsub equal substeps between the following points: 

Xi.^t = ^0 + (xi - xo)^^;^, isub = .. Nsub- (A-9) 

^^sub 

Then use ()A-8P to obtain C(xj^^^, k) and r]{xi^^^,k) from ^^{xi^^^-i, k) and r/(xi^^^_i, k). Here 
isub is the number of the substep. When xi is reached, remember the values ^(xi, k), r]{xi, k) 



X 

j exp [-A~'^Bs) ds 



A-^c. 



(A- 
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and increase Ng^b twice and do another iteration. Eventually, stop iterating on 'level-2'' 
after Nsub becomes large enough so that the resulting pair ^(xi, k), ri{xi, k) obtained at the 
current iteration is close enough to that from the previous iteration. 

What values should the integrator use at the Hevel-Z" iteration in the vector c to 
obtain y{xi^^^, k) from fc)? The easiest way would be to form c from ^{xi^^^^-i, k) 

and ri{xi^^^-i,k). But expecting this explicit method to have little stability, as typical 
of such methods, I decided to use an implicit method and to form c from (^(xi^^^,k) and 
r](xi_^^^,k). Of course, the code does not know the values at Xi^^^^ when it integrates from 
to which is what the explicit methods are all about. So I use a iterative 
procedure for that matter. First, assume that values of ^ and tj at the end of the substep 
are the same as at the beginning, and obtain the preliminary values at Xi^^^. These values 
are then used to repeat the substep as many times as it takes to get y{xi^^^) at the current 
iteration close enough to the one in the previous iteration. 

Throughout the solution, for I assume the following: 

• Quantities f«;(x), Pcr(2;,p), Tr{x), which are defined at x-grid planes, are inter- 
polated linearly in between the planes; 

• Spatial derivatives of the above quantities, u'{x), v'^{x), P^j.{x,p), are uniform between 
the grid planes. Their values correspond to the slopes of the linear interpolation of 
the above quantities; 

• Wave speed vq according to Equation (j4.3p : 

• The resonant wavenumber feres is related to the momentum pres as /cres^^ = 1- 
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Appendix B 

Diffusive flow incident on a moving absorbing boundary 



ulation as if they are crossing the position, at which they are placed, for the first time 
in their histories. The angular distribution of these particles is thus equal to the angular 
distribution of particles diffusively moving with respect to a flowing background medium, 
and incident on a fully absorbing boundary (the flow is directed into the boundary in our 
one-dimensional case). If the speed of the flow, u, is greater than the speed of the parti- 
cles, with respect to the flow, v, then, assuming isotropic distribution of particles in the 
reference frame tied to the flow, their angular distribution may be written as (j3.30p . This 
is simple, because for v < u all particles cross every position in the flow just once, because 
they cannot move upstream in the stationary reference frame. However, for v > u, this 
situation becomes more complicated, because particles are able to move forward as well as 
backward (against the flow) and the flux of such particles on a fully absorbing boundary is 
more difficult to estimate. 



upstream and propagating them downstream till they cross the fully absorbing boundary 
at X = for the first time. I recorded the angular distribution of these particles and fitted 
them with a simple scaling. I chose to search for the distribution in the power law form 



where v^f^ x is the x-component of the incident particle velocity measured in the shock frame 
(i.e., in the frame in which the absorbing boundary is at rest), fmin = 0, fmax = u + v, and 
C is found from condition 



It was mentioned in Section [3.1.51 that particles must be introduced into the sim- 



To solve this problem, I ran the Monte Carlo simulation, injecting the particles far 
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C 



a + 1 



(B-3) 
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making 



a + l 



— V, 




if Vj 



'mm 



'max; 



(B-4) 



'max 



mm 







otherwise. 



In order to derive a, one needs to minimize the following function of a to find the least 
squares fit: 



where the index i runs over all numerical bins of the speed Vgf^ x) the values Vi are the centers 
of these bins, and fi is the properly normalized fraction of incident particles that had the 
x-component of velocity in the i-th bin upon their incidence. I used a bracketing method 
to find the minimum, searching for it in the region a G [0.5, 2.0]. 



mono-energetic in the plasma frame, with a speed v, into a flow with speed u. The angular 
distribution that I used for these particles did not matter, because they were given enough 
time to scatter in the flow an isotropize before they reached the absorbing boundary. I 
covered the range v/u £ [1 ... 15]. For each such run, I fitted the angular distribution of 
particles first entering the shock with a single parameter power law and derived an a for 
this run. Then I plotted the resulting power law index a versus the ratio v/u of a run. The 
resulting curve a{v/u) can be described by the following simple equation: 



The angular distribution function (|B-4p with a given by ()B-6P is simulated in the 
code in order to introduce particles with a plasma frame speed v greater than the local flow 
speed u. Note that for v ^ u, the power law index approaches a — > 1.5, and for v ^ u the 
power law index approaches a — > 1.0, and it stays a = 1 for v < u, where (I3.30p becomes 
applicable. The last statement was demonstrated separately, in other simulations, and is 
obvious: for small v there are no backward-moving particles in the shock frame, so every 
particle crossing a plane crosses it for the first and the last time. 





In order to collect the data, I ran 30 simulations, introducing particles that were 





